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Abstract We discuss the physics and computation of lattice QCD, a space- 
time lattice formulation of quantum chromodynamics, and Kenneth Wilson’s 
seminal role in its development. We start with the fundamental issue of con¬ 
finement of quarks in the theory of the strong interactions, and discuss how 
lattice QCD provides a framework for understanding this phenomenon. A 
conceptual issue with lattice QCD is a conflict of space-time lattice with 
chiral symmetry of quarks. We discuss how this problem is resolved. Since 
lattice QCD is a non-linear quantum dynamical system with infinite degrees 
of freedom, quantities which are analytically calculable are limited. On the 
other hand, it provides an ideal case of massively parallel numerical compu¬ 
tations. We review the long and distinguished history of parallel-architecture 
supercomputers designed and built for lattice QCD. We discuss algorithmic 
developments, in particular the difficulties posed by the fermionic nature of 
quarks, and their resolution. The triad of efforts toward better understand¬ 
ing of physics, better algorithms, and more powerful supercomputers have 
produced major breakthroughs in our understanding of the strong interac¬ 
tions. We review the salient results of this effort in understanding the hadron 
spectrum, the Cabibbo-Kobayashi-Maskawa matrix elements and CP viola¬ 
tion, and quark-gluon plasma at high temperatures. We conclude with a brief 
summary and a future perspective. 
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1 Introduction 

In early 1974, Kenneth Wilson circulated a preprint entitled ” Confinement of 
Quarks”. The paper was received by Physical Review D on 12 June 1974, and 
was published in the 15 October issue of that journal pQ. In this paper, Wilson 
formulated gauge theories on a space-time lattic^] Using an expansion in 
inverse powers of the bare gauge coupling constant, Wilson demonstrated 
that lattice gauge theories at strong coupling confine charged states. He also 
argued that the absence of Lorentz invariance (or Euclidean invariance for 
the imaginary time used for the lattice formulation) is not a hindrance if 
there is a second order phase transition at some value of the gauge coupling 
constant, for in the vicinity of such a phase transition the lattice spacing can 
be taken to zero while fixing the physical correlation length at a finite value. 

This paper laid the conceptual foundation for understanding the quark 
confinement phenomenon. It showed that large quantum fluctuations of gauge 
fields at strong coupling can generate a force between charged states which 
stays constant at arbitrary distances. This is a novel type of force, essentially 
different from the Yukawa force due to exchange of particles which tends to 
zero at large distances. 

Initially, Wilson’s idea did not catch on rapidly. Techniques were hard 
to come by which allowed calculations of physical quantities such as hadron 
masses and connect them to physical predictions in the continuum space- 
time. The situation changed dramatically around 1979-1980 when Creutz, 
Jacobs and Rebbi [4], and Wilson himself 0, showed the possibility of nu¬ 
merically calculating the observables on a computer. Particularly dramatic 
was the calculation of the static quark-antiquark potential by Creutz for 
SU(2) gauge group in 1979 [6], and the calculations in 1981 of hadron masses 
by Weingarten (7j and by Hamber and Parisi ||. 

Traditional quantum field theory until that time could only deal with 
weakly coupled bound states such as hydrogen. The possibility of a technique 
which enables calculation of the properties of relativistic and strongly coupled 
bound states such as pion and proton was entirely new. 

The timing was also perfect from computational point of view. The CRAY- 
1 supercomputer which appeared in 1976 had revolutionized scientific com¬ 
puting, and lattice QCD computation could quickly exploit vector supercom¬ 
puters in the 80’s. Perhaps more important in retrospect, rapid development 
of microprocessors in the 70’s stimulated more than a few groups of particle 
theorists around the world to start developing parallel computers for lattice 
QCD. 

The development of lattice QCD has been continuous since then. With 
large-scale numerical simulations on parallel supercomputers, understanding 
of physics of lattice QCD progressed, which in turn led to better algorithms 
for computation. These algorithms allowed a better exploitation of the next 


1 Contemporary research toward lattice gauge theory was described by Wilson 
in his plenary talk at 1983 Lepton Photon Symposium [2], and in more detail in 
his historical talk at 2004 International Symposium on Lattice Field Theory at 
FNAL [3]. 
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generation of more powerful computers, which brought even more progress 
of physics. 

In the four decades of progress, lattice QCD has brought a deep under¬ 
standing on the physics of the strong interactions. It has matured to the point 
where one can make calculations with the physical values of quark masses, 
on lattices with sufficiently large sizes, at lattice spacings small enough so 
that a continuum limit can be carried out with confidence. 

In this article we review lattice QCD from four perspectives in the follow¬ 
ing four chapters. In Chapter |2j we discuss the foundation of lattice QCD, 
touching upon the connection between quantum field theory in Euclidean 
space time and statistical mechanics which was consciously exploited by Wil¬ 
son. We show how it led to a conceptual breakthrough in the understanding 
of confinement. We also explain the issues related with chiral symmetry. In 
Chapter[3j we discuss the computational aspects. We review how lattice QCD 
embodies an ideal case of massively parallel computation, and how this led to 
the development of parallel supercomputers for lattice QCD, which impacted 
seriously on the history of supercomputers up to the present time. Also dis¬ 
cussed is a special computational difficulty posed by the fermionic nature 
of quark fields, and how overcoming that difficulty has led to the algorithm 
in use today. In Chapter [4] we discuss some major physics results achieved 
so far in lattice QCD. The themes include the hadron mass spectrum, the 
determination of the Cabibbo-Kobayashi-Maskawa matrix and CP violation 
in the Standard Model, and the properties of quark gluon plasma at high 
temperatures and densities. Finally, in Chapter [5] we collect some thoughts 
on how Kenneth Wilson’s thinking and vision helped develop the subject. 

The tone of this article is partly historical, describing the development of 
lattice QCD and the role Kenneth Wilson played in it. It also reviews the 
achievements made in the four decades since its inception in 1974. 


2 Quantum chromodynamics on a space-time lattice 

2.1 Hadrons, Quarks, Quantum Chromodynamics 

If one looks up the Reviews of Particle Physics web page [9], one finds there 
an entire list of particles and their properties experimentally discovered to 
date. In addition to “Gauge and Higgs Bosons”, “Leptons” and “Quarks”, 
there are two lists named “Mesons” and “Baryons”, each of which contains 
hundreds of particles. Protons and neutrons, which make up atomic nuclei, 
are two representative particles belonging to the family of baryons. Pions and 
K mesons are less familiar, but important particles for binding protons and 
neutrons into nuclei, and they belong to the family of mesons. The mesons 
and baryons are collectively called “hadrons”. Their chief characteristic is 
that they participate in the strong interactions in addition to the electro¬ 
magnetic and weak interactions, while leptons participate only in the latter 
two interactions. 

Many of hadrons were discovered in the accelerator experiments in the 
50’s and 60’s. In 1964 Gell-Mann and Zweig proposed that hadrons are 
composed of more fundamental particles, which Gell-Mann named quarks. 
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Quarks were predicted to have unusual properties such as fractional charge 
in units of electron charge. Evidence has gradually built up, however, that 
quarks are real entities. Yet experimental efforts for detecting them in iso¬ 
lation have been unsuccessful. This situation is often called “quark confine¬ 
ment” . 

Quantum chromodynamics (QCD) proposes to explain the constitution 
of hadrons from quarks and their interactions. It is a quantum field theory 
with local SU(3) gauge invariance in which the quark field q(x), with x the 
space-time coordinate, transforms under the fundamental 3 representation of 
SU(3). The SU(3) quantum numbers are called color since the basic premise 
of the theory is that only color neutral states, trivial under SU(3), carry 
finite energy and hence exist as physical states. This is a general statement 
of ” quark confinement”. 

Local gauge invariance is a requirement that the frame of reference of an 
internal symmetry may be freely rotated at each space-time point without 
altering the content of the theory. Thus QCD as a gauge theory is to remain 
invariant under the transformation q(x) —?► V(x)q(x ) where V(x) E SU(3) 
may vary from point to point in space-time. This invariance requires the 
existence of a vector gluon field A fl (x) with values in the Lie algebra of 
SU(3) which tells how the local frame at a point x and at a different point y 
are related. 

QCD as field theory thus contains gluon as well as quark fields. It is 
defined by the Lagrangian density given by 

£qCD = ^2 Tr ( F ^( X ) 2 ) + v( x ) (n^ + m ) l( X )- (!) 

Here, q(x) = (qi(x)j • * •, qN f ) is a spin 1/2 Dirac spinor field for quarks, 
Nf denotes the number of quark flavors with m = (mi, • • •, rriN f ) the quark 
mass matrix, and D^ — iA^(x) is the covariant derivative with = 

d^A u — djyA^ + [A^, A v \ the gluon field strength, and g is the QCD coupling 
constant. 

The discovery of asymptotic freedom mm in 1974 showed that the 
coupling strength of non-Abelian gauge theories decreases toward zero at 
large momenta. Since deep inelastic electron nucleon scattering experiments 
carried out in the late 60’s indicated just such a behavior called scaling, 
the discovery boosted QCD to the leading candidate of the theory of strong 
interactions. 

A beautiful prediction of asymptotic freedom is the existence of logarith¬ 
mic violation of scaling which can be quantitatively calculated via renormal¬ 
ization group methods. The prediction was later confirmed by experiments, 
thus establishing the validity of QCD beyond doubt at high energies. 

Since asymptotic freedom at high energies means that the coupling strength 
increases in the opposite limit of low energies, it was natural to speculate that 
QCD also provided a solution to the long standing puzzle that quarks had 
never been observed in experiments. However, quantum field theory at the 
time, though quite sophisticated, did not possess means to analyze the be¬ 
havior of QCD for large coupling constant expected at low energies. 
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Space-time continuum 


Space-time lattice 



Fig. 1 Replacing continuum space-time by a space-time lattice. Quark fields are 
placed at lattice sites, and gluon fields on links. 


2.2 Formulation of QCD on a space-time lattice 

An essential ingredient for a strong coupling analysis is a mathematically 
well-defined formulation of QCD in which ultraviolet divergences are con¬ 
trolled. Kenneth Wilson approached this problem with three key ideas: (i) 
use Euclidean space-time with imaginary time rather than Minkowski space- 
time with physical time, (ii) replace Euclidean space-time continuum by a 
discrete 4-dimensional lattice to control ultraviolet divergence, and (iii) main¬ 
tain local gauge invariance as the guiding principle to construct field variables 
and action on the lattice. 

The use of Euclidean space-time brings out a beautiful and powerful con¬ 
nection between quantum field theory and statistical mechanics. This connec¬ 
tion was realized in the 60’s by Kurt Symanzik and others from the viewpoint 
of rigorously defining quantum field theory m- The explosive development 
in the theory of critical phenomena due to Leo Kadanoff, Michael Fisher, 
Wilson himself and others in the late 60’s and early 70’s brought the impor¬ 
tance of the concept to the foreground m • a formal proof that quantum 
field theory in real time can be recovered from that in imaginary time under 
a set of axioms was given in the first half of 70 ’sin. This connection, then, 
was a hot topic at the time and Wilson consciously exploited it in his work 
on lattice QCD. 

Let us consider a simple cubic 4-dimensional lattice in Euclidean space- 
time. The lattice points, called sites, are labeled by an integer component 4- 
vector n = (ni, n 2 , 713 , 714 ) G Z 4 , and have a physical coordinate x = na with 
a the lattice spacing. A pair of neighboring sites at n and n+jl with (1 the unit 
vector in the direction fi = 1, 2, 3,4 are connected by the link between the two 
sites which may be denoted as i = (n, /i), n G Z 4 , fi = 1, 2, 3,4. An elementary 
square, or plaquette, on the lattice is then labelled by P = (n, /i, u) with n 
the lowest corner and the two directions spanning the square denoted by /i 
and v. The lattice construction is illustrated in Fig. [lj 

It is natural to place a quark field q[x) on each site q[x) = q n with x = na. 
If one replaces the derivative d^q(x) by a finite difference ( q n +ii — q n -fi )/2n, 
the quark Lagrangian on a lattice will contain a bilocal term q n q n ±fi • 
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In the continuum space-time, bilocal quantities such as q(y)q(x) are ren¬ 
dered gauge invariant by inserting a path ordered phase factor U(y,x) = 
Pexpfy A tJb (x)dx fJL ), which transforms as U(y,x ) —)> V^(y)U(y, x)V{x) un¬ 
der the local gauge transformation. Wilson proposed to employ the phase 
factor along the lattice link U n/l ~ Pexp(i A /Ll (x)dx JLL ) as the funda¬ 

mental gluon field variable on the lattice. The lattice quark action can then 
be made invariant by inserting appropriate U n ^s to the bilocal terms. One 
thus finds for the quark action on the lattice, 

Squark = ^ \ Qnl^nfj,q n -\-fi ~ Q.nljl ^Qn-fij + a g n mpg n , (2) 

n/i n 

with mo the bare quark mass matrix. 

A classical action for gluons which reduces to the continuum action in 
the limit a —» 0 can be constructed by taking the product of U n/Jj : s around 
the boundary of a plaquette P [T]: 


-^gluori 


f'O p 


n u w > 

n/i£dP 


( 3 ) 


where Tr stands for trace over SU(3) indices. Here go is the bare gauge 
coupling constant at the energy scale of the lattice cutoff 1/a. 

Lattice QCD as quantum theory is defined by the Feynman path inte¬ 
gral. If 0(U, q , q) is an operator corresponding to some physical quantity, the 
vacuum expectation value over the quantum average is given by 


(0(U,q,q)) 


1 

Zqcd 




dq n dq n O(U, q, q) exp(S Q c D ) 


( 4 ) 


with 


Zqcd 




dq n dq n exp(5 QC D)- 


( 5 ) 


where 


^QCD — ^gluon T ^quark- (6) 

The gluon link variable U nM takes values in the group SU(3) rather than 
the Lie algebra for the vector gluon field The integration over U n/1 ’s 

should be defined as invariant integration over the group SU(3). Since SU(3) 
has a finite volume under this integration, the lattice Feynman path integral 
is well-defined without gauge fixing. 

The integration over the quark fields also requires some care. Since quarks 
are fermions, the path integral has to be defined in terms of Grassmann 
numbers which anticommute under exchange, i.e., q n qm — — qmqn • These 
points are clearly spelled out in the Wilson’s original paper in 1974 [I] . 
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R 


R 



Fig. 2 Left panel: Wilson loop for a temporal size T and a spatial size R measures 
the potential energy for a pair of static quark and antiquark. Right panel: The 
leading order of the strong coupling expansion is obtained by tiling the Wilson 
loop by plaquettes. 


2.3 Confinement of color 

Whether quarks are confined in QCD can be examined if one knows how 
to calculate the energy of an isolated quark in interaction with the gluon 
fields: an isolated quark would exist only if the energy of that state is finite. 
An elegant method invented by Wilson is to consider a pair of static quark 
and antiquark, which is created at some point in space-time, then separated 
to some distance, stays in that configuration for some time, and finally is 
brought together to a point and annihilated. Geometrically, the space-time 
trajectory of the pair forms an oriented closed loop C. Since the color charge 
of the quark interacts with gluon fields, the creation and annihilation of the 
pair inserts a phase factor 

W [< C] = P exp j> dx^A^x^j (7) 

in the path integral where P indicates ordering along the path. This is the 
Wilson loop operator. 

If the loop has the shape of a rectangle of a width R in the spatial direction 
and an extent T in the temporal direction as shown in the left panel of Fig.[2j 
the quantum average of the Wilson loop operator measures the energy E(R) 
of the quak-antiquark pair relative to the vacuum over a temporal length T 
so that 


(W [R x T]) oc exp (-E(R)T ), T oc. (8) 

Wilson speculated that for large values of the gauge coupling, fluctuations 
of gluon fields would be large, leading to a significant cancellation among the 
contributions to the Wilson loop. This would mean that the Wilson loop 
average rapidly decreases for larger loops, or equivalently the probability 
decreases that the quark and antiquark are found in a well-separated config¬ 
uration. This would mean confinement. Another equivalent statement would 























be that the energy E(R) grows for larger separations R , so that separating 
quark from antiquark is not possible. 

An amazingly simple calculation suffices to verify this picture if one em¬ 
ploys lattice QCD. Let us assume, following Wilson, that quark fields do not 
play an essential role. When the bare coupling go becomes large, the gluon 
path integral in @ can be calculated by an expansion in inverse powers of 
g$. The leading term is obtained when one tiles the R x T surface of the 
Wilson loop by a set of plaquettes from the expansion of the gluon part of 
the weight exp (S g i uon ) : as shown in the right panel of Fig [ 2 ] We then find 
that 

/ 1 \RT/a 2 

(W[RxT]) = N c ^- w -iJ (1 + OQ? 0 - 8 )), gl ^ oo, (9) 
with N c = 3 for SU(3) so that 

E(R) -» aR, R-> oo with a = L (log (N c gl) + 0(g^ 8 )) , (10) 

namely, a static pair of quark and antiquark are bound by a potential linearly 
rising with the separation. Hence they cannot be separated to infinite distance 
with any finite amount of energy. 

A closed loop C has two geometrical characteristics, the length of the 
loop L[C] and the area of the minimal surface A[C] spanned by the loop. 
The confinement behavior corresponds to an area decay for large loop; 

(W [C]} oc exp (— aA[C]) . (11) 

If, on the other hand, the Wilson loop expectation value decays with the loop 
length, 


(W [C]) oc exp (—/iL[C ]), 


( 12 ) 


the energy of a quark-antiquark pair saturates to a constant E(R) g 
for large separation R oo. Hence there is no longer confinement. The 
confining and non-confining phases of gauge theories are thus distinguished 
by the behavior of the Wilson loop expectation value. 

The possibility that there can be both confining and non-confining phases 
raises an interesting question whether a confining phase can turn into a non¬ 
confining phase if some parameters of the theory are varied. Temperature 


is such an important parameter. As we discuss in more detail in Sec. 4.4 


the confining property becomes lost through a phase transition when the 
temperature is raised sufficiently. 


2.4 Continuum limit 


Let us go back to the strong coupling result in (10). If one uses the phe¬ 
nomenological value y/a ~ 440MeV, and if one assumes that a value log (7V c $o) 
0(1) is sufficient for the strong coupling expansion to converge, the corre¬ 
sponding value of the lattice spacing equals a ~ 0.5 fm with 1 fm= 10“ 15 m. 
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Fig. 3 Monte Carlo result for the static potential energy E(R) calculated in pure 
gluon theory showing a linearly rising term at long distances and a Coulomb term 
at short distances. From Ref. p~5] . 


This value is comparable to a typical length scale of the strong interactions, 
e.g., the charge radius of proton R p ~ 0.9 fm. We certainly know that space- 
time is continuous well below such length scales. Thus one need to know if 
confinement holds for smaller lattice spacings. 

In Fig. [3] we show the static potential E(R) calculated by Monte Carlo 
simulation of the pure gluon theory |T5j. We observe a linearly rising potential 
E(R) (jR at large distances. There is also a Coulomb behavior V(R) ~ a/R 
at short distances, which is consistent with a perturbative one gluon exchange 
expected from asymptotic freedom. The lattice spacing estimated from the 
phenomenological value y/a ~ 440 MeV is a ~ 0.0544(4) fm. This is one of 
many evidences that the property of confinement holds not just at strong 
coupling but also toward weak couplings with small lattice spacings. 


The question then is whether the confinement property really persists 
as the lattice spacing is taken toward the limit of continuous space-time 
a —>> 0. For Wilson, who elucidated critical phenomena with his renormal¬ 
ization group ideas, this was not a conceptually difficult issue. The strong 
coupling calculation demonstrates that the system of gluons without quarks 
is in a confining phase. When one decreases the coupling constant go , this 
property persists as long as one does not encounter a phase transition to 
a non-confining phase characterized by perimeter decay (12). Suppose that 
there is such a second-order phase transition at g^ = g%. Close to the critical 
point, the correlation length measured in lattice units diverges £(go)/a 00 

as go g c • This means that one should be able to hold the physical cor¬ 
relation length £(go) fixed to a physical value observed in experiment while 
sending the lattice spacing to zero a —»• 0, recovering physics in the space-time 
continuum. 
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One can define lattice gauge theory for a variety of groups and space-time 
dimensions. As with the case for spin systems, the existence and position of 
second-order phase transition points would depend on them. For the group 
SU(3) in 4 dimensions, a particularly attractive possibility is g c = 0. In fact, 
this is the only possibility if confinement at low energies and asymptotic 
freedom at high energies are to coexist in the same phase. 

The renormalization group allows a more concrete discussion. The change 
of gauge coupling go —>• go+dgo under a change of cutoff scale a a-\-da while 
fixing a physical scale constant defines a renormalization group /? function 

a dg^a)_ = _p( go ^ (13) 


where the minus sign is inserted to be compatible with the conventional 
definition using momentum slicing. The strong coupling result (10) shows 
that /3(go) is large and negative for large go. For small values of go , one can 
employ perturbation theory to confirm the asymptotic freedom result, 


with 


P(9o) = “Mo ~ Mo + O(g 7 0 ) 


b 0 = 

h = 


1 


— N r - -N 


(47r) 2 ^ 3 

1 ( 34 at 2 Ac 2 - 1 10 Ar 

t — — — iVr --- N f - N c N f 

(47r 4 V 3 N c } 3 f 


(14) 

(15) 

(16) 


for SU(N C ) gauge group and Nf flavors of fermions in the fundamental rep¬ 
resentation. The first two coefficients are negative for our world with N c = 3 
and Nf = 6. A natural supposition is that the beta function is negative for 
all values of the coupling, with the only zero residing at g c = 0. 

Wilson started a numerical Monte Carlo calculation to check this by a 
block spin renormalization group for gauge group SU(2) [5]. This attempt 
was followed by several serious calculations with SU(3) gauge group in the 
80’s. The results, though mostly restricted to the pure gluon theory and 
had fairly large errors, showed that the beta function stayed negative and 
became consistent with the two-loop result (14) toward weak coupling (see, 
e.g., Refs. [MED- 

A different, and perhaps a more elegant, approach was developed in the 
90’s by the Alpha Collaboration m- Their method of step scaling func¬ 
tion defines the renormalized coupling constant at a scale 1/L through the 
Schrodinger functional on a lattice of a finite size L 3 x T with a prescribed 
boundary condition in time, and follows the evolution of the coupling under 
a change of scale by a factor 2. The continuum limit is systematically taken 
in this process. Therefore, the end result is the evolution of the renormalized 
coupling from low to high energies in the continuum theory. 

In Fig. [4] we show the result for the pure gluon theory (solid circles) 
and a comparison with perturbation theory (dashed and dotted lines). The 
full evolution runs from 10 GeV down to a few hundred MeV, where the 
system is in a strong coupling non-perturbative regime. Thus the confining 
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Fig. 4 Running of renormalized coupling constant ot(q 2 ) = g 2 (g 2 )/47r in the 
Schrodinger functional scheme as a function of momenta scale q. Dashed and dot¬ 
ted lines are two- and one-loop evolution starting from the right- most point. The 
physical scale is fixed by the Sommer scale ro = 0.5 fm. From Ref. [T 8 ] . 

behavior at large distance is continuously connected with the asymptotically 
free behavior at short distances. The full evolution runs almost parallel to 
the two-loop evolution. Similar results have been obtained in full QCD with 
Nf = 2,3,4 dynamical quarks [E312QEE]. 

Based on the results described above, one can state, though a mathemat¬ 
ically rigorous proof is yet lacking, that QCD is in the confining phase over 
the entire range of coupling from zero to infinity. 


2.5 Quarks and chiral symmetry 
2.5.1 Chiral symmetry 

An important feature of the strong interaction which was recognized in the 
late 50’s and early 60’s is chiral symmetry. Studies in this period led to the 
concept of spontaneous breakdown of symmetry, and a realization that it is 
accompanied by the emergence of massless bosons, called Nambu-Goldstone 
bosons today. 

In QCD language, chiral symmetry is invariance under the global trans¬ 
formation q(x) Wq(x ), where W = exp(ia^) £ SU(Nf) acts on the 

Dirac-flavor indices and rotates left and right handed chiral components of 
q{pc) in the opposite direction. This symmetry is explicitly broken by the 
quark mass term. Hence it holds approximately for three light quarks, up, 
down and strange. It is not very relevant for the heavier quarks, charm, bot¬ 
tom and top. The octet of pseudo scalar mesons, 7 r, K, 77 , are identified as the 
Nambu-Goldstone bosons corresponding to spontaneous breakdown of chiral 
symmetry. With a strongly interacting dynamics at large distances, QCD 
should dynamically explain spontaneously broken chiral symmetry and its 
physical consequences. 
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2.5.2 chiral symmetry on a space-time lattice 

Introduction of a space-time lattice affects bosons and fermions in different 
ways. This is most easily seen by looking at the kinetic term in the free field 
case. For a boson field </> n , the second-order derivative —d 2 (/)(x) is discretized 
as — {4>n+ji + ~ 2 </> n )/ 2 a 2 , whereas for a fermion field q{x ), the first or¬ 

der derivative 7 ^d^q^x) is discretized as 7 M (g n+ ^ —g n _^)/2a. The momentum 
space expression then becomes 2(1 — cos (p^a))/a 2 and 27 ^ sin(p M a), respec¬ 
tively. There is a crucial difference in the location of the zeros in the Brillouin 
zone: the boson case has a zero only at p^ = 0 for all /i whereas the fermion 
case has a set of zeros for p^ = 0 or ir/a for each fi. Since the zeros gives rise 
to poles in the propagator, a naive fermion discretization leads to multiple 
copies of the state at p M = 0. This is the fermion species doubling problem. 

Building upon a pioneering work by Karsten and Smit [ 22 ], Nielesen and 
Ninomiya m proved that essentially the same conclusion holds under a set 
of rather general conditions. The theorem states that if the fermion action 
satisfies (i) chiral symmetry q n —>> exp (icr /5 )g n , (ii) invariance under unit 
translation on the lattice, (iii) Hermiticity, and (iv) locality, then the spec¬ 
trum of fermions contains even number of particles, half of them left handed 
and the other half right handed. 

An elegant topological proof runs as follows [21]. If we write a Dirac 
fermion action on a lattice in a general form, Sf = the 

assumptions (i) and (ii) imply that D n;m = Yh^lnF^n — m) where F^(n) is 
a vector function which, by (iii) satisfies F^(n) = F’ /x (—n), and by (iv) rapidly 

decreases as \n\ becomes large. The Fourier transform F^(p) is, therefore, a 
well-defined real vector field over the Brillouin zone in momentum space. 
Let p = = 1, 2, • • • be the zeros of the vector field F^(jp) = f^ u {Pu — 

Pv^) + 0((p—p^) 2 ). They correspond to the poles of the propagator D _ 1 (p), 
and hence to the particle states. The relative chirality of these states are 
determined by the index, sign(det f^). Since the sum of index of a vector 
field over the 4-dimension torus vanishes by the Poincare-Hopf theorem, there 
has to be an equal number of fermion states with opposite chirality. 

The Nielsen-Ninomiya theorem indicates that one either has to abandon 
chiral symmetry or one has to allow for the presence of species doubling. 
Wilson’s choice, which he actually wrote down [25! prior to the publication 
of the theorem, was to add a term which softly breaks the chiral symmetry 
of the naive lattice action in ®): 

a 3 / \ 

S $W = y Yi ( 2 9n<?n - q n Un^qn+ii ~ ■ (17) 

nil 

For the free field case, this term adds JE (1 — cosp^a)/a to the kinetic term 
i 7 ii sinp^a, and hence removes the zeros at p^ = n/a. 

Let us add that toward the continuum limit a 0 the Wilson’s added 
term becomes of form aq(x)D^q(x) relative to the original term q{x)i^(^q{x) 
in ©, *.e., higher order in a. Apart from removing the doublers, the effect 
of the added term disappears in the continuum limit. 
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Chiral breaking effects of the Wilson’s added term can be analyzed by the 
method of Ward identities [26]. In particular, a definition of quark mass can 
be given that satisfies the PCAC relation nm. a detailed analysis of the 
phase structure on the (go, mo) plane was made [28], and the existence of a 
massless pion, in spite of chiral symmetry breaking, was explained as due to 
spontaneous breakdown of parity Z(2) symmetry [29] . With these analytical 
developments, Wilson’s formulation provides a quantitative framework for 
the computation of physical observables, and is used extensively in Monte 
Carlo studies. 

A variant of Wilson’s formulation is to consider two flavors of quarks as 
a pair and add a twisted mass term to the naive action © m : 


5S tm = a 4 '^2q n iVq'y5T3qn, 


(18) 


n 


where 73 acts on the flavor index. This formulation has an attractive feature 
that one can twist the angle a = arctan g q /mo in such a way that 0(a) lattice 
artifacts are absent in physical observables |BIUB2] . Large-scale simulations 
are being made using such ’’maximally” twisted QCD. 

There is a different method, called the staggered formulation [33], which 
retains a U(l) chiral symmetry at the cost of a four fold species doubling. 
In the 4-dimensional Euclidean formulation [34] , it starts with a single com¬ 
ponent fermion field Xn at each site, and reconstructs four species of Dirac 
fields if = (Vq, • • •, Vq) from 16 y’s on the 16 vertices of a 4-dimensional 
unit hypercube [35] • The original action is invariant under an even-odd U(l) 
symmetry Xn exp(i(—l)l n l<a)x n , which translates into an axial U(l) trans¬ 
formation on the four Dirac fields of form ^ —>> exp (ia^s 0 ^ 5 )^ where £5 acts 
on the species index a of the Dirac field Vkn a = 1, • • •, 4. 

It is generally believed that lattice QCD with the staggered fermion for¬ 
malism converges to continuum QCD with Nf = 4 degenerate flavors of 
quarks [36 . An elaborate effective theory description has been developed 
to control the breaking of the full 4 ’’flavor” chiral symmetry down to the 
U(l) subgroup at finite lattice spacing [37] . On these theoretical bases, the 
staggered formalism is also extensively used in Monte Carlo simulations. 

2.5.3 Lattice fermion action with chiral symmetry 

In 1981, a year after Nielsen and Ninomiya presented their theorem, Wil¬ 
son revisited the issue of chiral symmetry from a different perspective with 
Ginsparg [38] . He asked what would be the relation satisfied by a lattice Dirac 
operator D n?m if it were derived by a (chiral symmetry breaking) block spin 
transformation from a chiral invariant theory. The answer turned out to be 
remarkably simple; it is given by 



(19) 


k 


Since D does not anticommute with 75 , assumption (i) of the Nielsen-Ninomiya 
theorem is not satisfied, and so the conclusion of the theorem does not hold. 
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In fact, there is no species doubling. Furthermore, the axial vector current 
for this action has the correct U(l) anomaly. 


One can rewrite (19) in terms of the propagator G = D 1 as 

+5^71,m T Gn,m'y5 — ^^nm'T5* 


( 20 ) 


Hence the breaking of chiral symmetry is a local effect. One may expect then 
that a modified form of symmetry may exist. This was found almost 20 years 
later [39]. Fermion actions that satisfy the Ginzparg-Wilson relation (19) are 
invariant under an infinitesimal transformation given by 


Sip = ip{ 1 - laD) 75 , Sip = 75(1 - 1 aD)ip. 


( 21 ) 


Several forms of fermion action which satisfy the Ginsparg-Wilson relation 
were discovered in the 90’s. One form is a domain-wall formalism [401 141 ] in 
which the 4-dimensional fermion field is constructed as the zero mode of a 5- 
dimensional theory generated by a mass defect at the boundary in a fictitious 
fifth dimension. Another form is given by the overlap formalism [421143] . In 
this case an explicit form of the operator D is given by 

D= 1 (l + 7 5#(# t #)“ 1/2 ) , H = 75 (aDw - 1) (22) 

with aDw + rno the operator for the Wilson fermion action with a negative 
mass mo = —1. The two forms are equivalent in the limit of infinite fifth 
dimension @145]. 

Yet another form is the perfect action @ 6 ], so named because it is defined 
as the fixed point of a block spin transformation of renormalization group for 
QCD. This form follows from the line of reasoning of Ginsparg and Wilson, 
but it was pursued and arrived at independently almost two decades later. 

All forms of action, particularly the domain wall and overlap actions, have 
come to be used extensively in the last decade. The domain wall formalism 
has been exploited by RBC-UKQCD Collaborations (see, e.g., m ), and the 
overlap formalism by JLQCD (see @ 8 ] for a recent review). The situation 
with the perfect action is reported in [49] . 


2.5.4 Spontaneous breakdown of chiral symmetry 

Spontaneous breakdown of chiral symmetry is best studied by examining 
the behavior of the order parameter of chiral symmetry. This is given by the 
quark bilinear operator D = (q n q n )• If A 7 ^ 0 after sending the spatial volume 
to infinity V -4 00 followed by the limit of quark mass to zero m q —)> 0 , then 
chiral symmetry is spontaneously broken. 

In Fig. [5] we show the result for A as a function of the degenerate up-down 
quark mass m uc [ in lattice units obtained with the overlap formulation for 
Nf = 2 + 1 lattice QCD. We choose these data since the overlap formalism 
is the cleanest regarding the chiral aspect among lattice fermion formula¬ 
tions. The terminology Nf = 2 + 1 refers to the fact that the up and down 
quark masses are taken degenerate, while the strange quark mass has a sep¬ 
arate value. Strictly speaking, the values shown in Fig. [5] are obtained from 
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Fig. 5 Chiral order parameter E — {q n q n ) as a function of degenerate up-down 
quark mass in lattice units in Nf =2 + 1 flavor lattice QCD with overlap fermion 
formalism. The strange quark mass is held fixed at a value close to experiment. 
From Ref. [48]. Color online. 


the eigenvalue distribution of the Dirac operator and not by calculating the 
condensate directly. 

The concave curvature as a function of m u d is consistent with the presence 
of a logarithm term predicted by chiral perturbation theory. Extrapolating 
to m u d = 0 including the effect of the logarithm yields a non-zero value, 
supporting spontaneous breakdown of chiral symmetry. The dependence of 
the pion mass is consistent with oc m u d (up to logarithmic corrections), 
as follows from the Nambu-Goldstone theorem. 

Similar results have been obtained for the staggered fermion formalism 
which has U(l) chiral symmetry. The analysis is more complicated for the 
Wilson fermion formulation, since one needs to carry out a subtractive renor¬ 
malization for both chiral condensate and quark mass due to the soft chiral 
symmetry breaking induced by the Wilson term (17). Once these subtrac¬ 
tions are done, one finds the signatures expected for spontaneous breakdown 
of chiral symmetry, such as the relation oc mn u d for an appropriately 
defined quark mass mm- 


2.6 Heavy quarks on a lattice 

The three light quark quantum numbers, z.e., up, down and strange, have 
been known from 1930’s and 1950’s. In contrast, heavier quarks were discov¬ 
ered more recently, i.e., charm quark in 1974, bottom quark in 1977, and 
top quark in 1995. These heavy quarks occupied the central place in the ex¬ 
perimental and theoretical studies toward establishing the Standard Model, 
particularly with the construction of B factories in the 90’s and experiments 
with them in the 2000’s. 

Studying heavy quarks with lattice QCD poses the problem that for a 
large value of a heavy quark mass rrih , the dimensionless combination ra^a 
is also large, leading to an amplification of lattice discretization errors, espe- 
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daily if m^a 1. This situation applies to the bottom quark with a mass 
mb & 5 GeV, since lattice QCD simulations to date have been made with 
inverse lattice spacings in the range a -1 ^2 — 4 GeV. 

Several methods have been formulated and employed to deal with this 
problem. The static approximation [50] is an expansion in l/m^, NRQCD [51] 
is a reformulation of QCD for non-relativistic quarks with an expansion in 
powers of the quark velocity u, and a relativistic formalism [5211531154] modifies 
the Wilson quark action so as to systematically reduce the effects of large 
m h a. 

All three methods have been extensively used to calculate physical quan¬ 
tities involving charm or bottom quark. In particular, matrix elements such 
as the pseudoscalar decay constants and form factors calculated with these 
methods have been playing an important role in constraining the Cabibbo- 
Kobayashi-Maskawa matrix elements including the CP violation phase. 

It should be mentioned that with increasingly smaller lattice spacings 
becoming accessible with progress of algorithms and computer power, direct 
simulations are replacing calculations with effective heavy quark theories. 
This has already occurred for charm quark, and it may not be too far into 
the future that bottom quark becomes treated in a similar way. 


3 Lattice QCD as computation 

3.1 Numerical simulation and lattice QCD 

Lattice QCD offered a framework for conceptually understanding the dynam¬ 
ics of non-Abelian gauge fields. In particular, it elucidated the mechanism of 
confinement in a way which would have been impossible in the perturbative 
framework of field theory. Nonetheless, calculation methods to obtain phys¬ 
ical results did not really exist; for example, higher order strong coupling 
expansions were very cumbersome and hard to extrapolate to the continuum 
limit expected at go 0. 

Monte Carlo simulations offered a new approach to solve this impasse. 
Of course Monte Carlo methods had been known since the pioneering era of 
electronic computers in the late 40’ and early 50’s. The Metropolis algorithm 
to handle multi-dimensional integrations for statistical mechanical systems 
was formulated in 1953 [55]. Applications to spin systems in statistical me¬ 
chanics started to appear in the 60’s and was pursued increasingly in the 
70’s. 

It is in this context that Creutz, Jacobs and Rebbi carried out a Monte 
Carlo study of Ising gauge theory with Z(2) gauge group in 4 dimensions in 
1979, finding a first-order phase transition [4]. Creutz extended the applica¬ 
tion to SU(2) gauge group |6., and extracted the string tension a in front 
of the area decay of the Wilson loop. The dependence of a on the gauge 
coupling go turned out to be consistent with the scaling law predicted by 
the renormalization group. This suggested the possibility that a continuum 
limit could be successfully taken, giving rise to a hope that the confinement 
problem could be solved in a numerical way. 
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Wilson’s interest in numerical analyses and computer applications started 
early in his career [3]. He often used numerical methods to carry out his 
analyses 0 In his 1979 lecture at the Cargese Summer School [5], Wilson 
reported a block-spin renormalization group analysis for SU( 2 ) gauge group 
using Monte Carlo methods to evaluate the necessary path integral averages. 

These attempts introduced a method hitherto unknown in field theory. 
The method looked very promising, and immediately attracted attention of 
particle theorists. In particular, it was applied to calculate masses of hadrons 
directly from quarks and gluons. The method was based on the observation 
that if On(t) is an operator at a time slice t corresponding to a hadron H, 
e.g., 7r + ~ u'jsd for pion or p ~ {^uC^d)d for proton, then the 2 -point 

Green’s function for this operator behaves for large times as 

^0h(£)Oh(O)^ ~^ Z ex P (—muat) , t oo, (23) 

where ran denotes the mass of hadron H. Therefore, calculating the two- 
point function by a Monte Carlo simulation and extracting the slope of the 
exponential decay would yield the mass of hadron H. 

This program was carried out by Weingarten [7], and by Hamber and 
Parisi [ 8 ] in 1981. The lattice size employed was 12 4 for the Icosahedral 
subgroup of SU( 2 ) for the former, and 6 2 3 x 12 for SU(3) for the latter work. 
Albeit very modest in today’s standards, their studies clearly demonstrated 
the feasibility of their approach, which accelerated an explosive development 
of Monte Carlo simulations in lattice QCD. 


3.2 Massive parallelism and lattice QCD 

Monte Carlo simulation of lattice QCD is based on the method of importance 
sampling. Let 0 n be a field variable at a site n on a lattice A with V lattice 
points, and 0(0) an operator. In field theory, one wishes to calculate the 
average value of 0 ( 0 ) weighted with the action 5[0] according to, 

(om = 4 / II exp (S[4>]), (24) 

** n 


where 


z=/n^°w ex pwi)' ( 2s ) 

^ n 

Let us define a configuration C as a point C = (0i, • • •, 0y) in the in¬ 
tegration space of all fields on the lattice A. Monte Carlo evaluation of the 

2 Probably his most fam ous numerical work is a renormalization group solu¬ 
tion of the Kondo problem [56] . The numerical rigor he maintained for this work 
has become legendary. For the universal ratio of the two temperatures Tk and 
To characterizing the high and low temperature scales, he obtained W/(47 t) = 

(47r) _1 Tic/To = 0.1032(5). Six years later, an exact solution by the Bethe ansatz 
yielded W/{ 47r) — 0.102676- • • 57], verifying the Wilson’s number to the fourth 
digit within the estimated error! 
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Space-time lattice 



Processor array 



Fig. 6 Mapping of space-time lattice A = U^^Aa to an array of processors 
P a ,<a = interconnected according to the space-time lattice topology. 

Color online. 


integral proceeds by setting up a stochastic chain of configurations C% 

C 2 Cn such that the distribution of configurations converges to the 

weight of the integral: 

Pn(C) = Jf E S (° - C 0 ->■ 6XP ( | [C]) » N ^°°- (26) 

i =1 

In the well-known Metropolis algorithm, the convergence is guaranteed by 
accepting or rejecting a new trial configuration C tr i a i according to the prob¬ 
ability P acc = Min( 1 ,exp(P[Ctri a i] — S'fC'oid])) where C Q id is the latest con¬ 
figuration of the chain. 

A very important property in Monte Carlo calculations of field theoretical 
systems is locality. A field variable <fi n at a site n interacts only with those 
variables </> m ’s in a limited neighborhood of the site n. In creating a trial 
configuration C tr i a i from an old configuration C Q id, and in calculating the 
difference of the action 5[Ct r i a i] — S'[GW], the numerical operations at a site 
n can be carried out independently of those at a site m unless the pair is 
within the limited neighborhood of each other. 

As shown in Fig. [ 6 j let us divide the entire lattice A into a set of sub¬ 
lattices A a , a = 1 , 2, • • •, Np, and assign each sub-lattice A a to a processor 
P a . Because of locality, calculations by the processor P a can be carried out 
independently of those by the other processors, except that the processors 
with overlapping boundaries have to exchange values of 0 * 2 s in the boundaries 
before and/or after the calculations in each sub lattice. This means that, for a 
fixed lattice size, the computation time can be reduced by a factor Np , and for 
a fixed sub-lattice size, one can enlarge the total lattice size proportionately 
to the number of processors Np without increasing the computation time. 

This is an ideal case of the parallel computation paradigm. The locality 
property, which is one of the fundamental premises of field theoretic descrip¬ 
tion of our universe, allows a mapping of calculations on a space-time lattice 
to a parallel array of processors interconnected with each other according to 
the connection of space-time sub-lattices. 
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3.3 Parallel computers for lattice QCD 

Immediately after Monte Carlo calculations started in lattice QCD, several 
groups started to plan the building of a parallel computer for lattice QCD 
calculations. A crucial factor which helped push such an activity was a rapid 
development of microprocessors in the 70’s. As shown in Table [l] starting 
with a 4 bit Intel 4004 in 1971, increasingly more powerful microprocessors 
were developed and were put out on the market at prices affordable by aca¬ 
demic scientific projects. 

In Table [2] we list the parallel computers developed by physicists for lattice 
QCD in the 80’s. Typically, these machines employed a commercial micropro¬ 
cessor such as those listed in Table [l] as the control processor, and combined 
it with a floating point unit to enhance numerical computation capabilities. 

The famous CRAY-1 vector supercomputer came on the market in 1976. 
Vector supercomputers developed rapidly, and dominated the market in the 
70’ and 80’s. However, the progress of parallel computers was even faster. By 
the end of the 80’s, parallel computers caught up and even overtook vector 
computers in speed. We can observe this trend by tracking blue and green 
symbols (vector computers) and red and violet symbols (parallel computers) 
from the 80’s to early 90’s in Fig. [7| 


Table 1 Representative microprocessors developed in the 70’s. Price/chip is from 
Wikipedia. 


year 

name 

bit 

price/chip 

1971 

Intel 4004 

4 bit 

$ 60 

1972 

Intel 8008 

8 bit 

$120 

1974 

Intel 8080 

8 bit 

$360 

1974 

Motorola 6800 

8 bit 

$360 

1978 

Intel 8086 

16 bit 

$320 

1979 

Motorola 68000 

32 bit 



Table 2 Parallel computers dedicated to lattice QCD developed in the 80’s. ’’year” 
marks the completion date. 


name 

year 

authors 

CPU 

FPU 

peak speed 

PAX-32* 

1980 

Hoshino-Kawai 

M6800 

AM9511 

0.5 MFlops 

Columbia 

1984 

Christ-Terrano 

PDP11 

TRW 

- 

Columbia-16 

1985 

Christ et al 

Intel 80286 

TRW 

0.25 GFlops 

APE1 

1988 

Cabibbo-Parisi 

3081/E 

Weitek 

1 GFlops 

Columbia-64 

1987 

Christ et al 

Intel 80286 

Weitek 

1 GFlops 

Columbia-256 

1989 

Christ et al 

M68020 

Weitek 

16 GFlops 

ACPMAPS 

1991 

Mackenzie et al 

micro VAX 

Weitek 

5 Gflops 

QCDPAX 

1991 

Iwasaki-Hoshino 

M68020 

LSI-logic 

14 GFlops 

GF11 

1992 

Weingarten 

PC/AT 

Weitek 

11 GFlops 


*not for lattice QCD 
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Table 3 Parallel computers developed for lattice QCD in the 90’s and later, ’’year” 
marks the completion date. 


name 

year 

authors 

CPU 

vendor 

peak speed 

APE100 

1994 

APE Collab. 

custom 

- 

O.lTflops 

CP-PACS 

1996 

Iwasaki et al 

PA-RISC 

Hitachi(SR2201) 

0.6TFlops 

QCDSP 

1998 

Christ et al 

TI DSP 

- 

0.6TFlops 

APEmille 

2000 

APE Collab. 

custom 

- 

0.8Tflops 

QCDOC 

2005 

Christ et al 

PPC-based 

IBM(BG/L) 

lOTFlops 

PACS-CS 

2006 

Ukawa et al 

Intel Xeon 

Hitachi 

14TFlops 

QCDCQ 

2011 

Christ et al 

PPC-based 

IBM(BG/Q) 

500TFlops 

QPACE 

2012 

Wettig et al 

PowerXCell 

- 

200TFlops 



Year 


JPN project machines 

• NWT 

• CP-PACS 

■ Earth Simulator 

• K 

Vector machines 
o CRAY/CDC 

• Hitachi/Fujitsu/NEC 
Vector parallel machines 

• Fujitsu 

■ NEC 
o CRAY 

Parallel machines 
o CRAY 

□ IBM 
O Intel 
A TMC 
☆ China 

• Fujitsu 

• Hitschi 

■ NEC 

QCD parallel machines 

□ Tsukuba 
ffi Columbia 
a APE 

N GF11 (IBM) 


Fig. 7 Peak speed of supercomputers. The circle at the bottom left is CRAY-1. 
Red fancy symbols show parallel machines developed for QCD. Color online. 


In Table [3] we list the parallel computers developed in the 90’s and later for 
lattice QCD. Big success continued with CP-PACS in Japan, which occupied 
the top position in the Top 500 list of supercomputers in November 1996, and 
QCDSP in USA. We observe an increasing involvement of major vendors such 
as Hitachi and IBM. This was necessary to secure advanced technology and 
computer building knowhow to develop those fast computers. Probably the 
most well-known of this trend is the QCDOC project, which gave rise to the 
IBM BlueGene/L, and the QCDCQ project which ran parallel with the IBM 
BlueGene/Q development. In this way lattice QCD has seriously influenced 
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the development of parallel supercomputers for scientific computing. Another 
trend one can observe in Table [3] is the use of commodity processors such 
as Intel Xeon for a quick machine building. This was the option adopted by 
PACS-CS and QPACE. 

Today, the fastest supercomputers have reached the peak speed of 0(50) 
Pflops [58]. This has been achieved in two ways. The K computer in Japan 
and Sequoia (BlueGene/Q) in USA connected O(10 3 * 5 ) multi-core processors 
running at 0(100) Gflops. Tianhe-2 in China and Titan in USA boosted the 
computing power by adding multiples of GPGPU’s running at 0(1) Tflops 
to each node. Further increase of computing speed faces a serious issue that 
the memory cannot supply data fast enough to the processing units and the 
power consumption is becoming too large for large systems. Serious effort is 
already going on, however, to overcome these problems. 


3.4 Fermion problem and hybrid Monte Carlo method 


Monte Carlo calculation for the gluon fields, though somewhat complicated 
by the SU(3) nature of the field variable, is straightforward. Calculations for 
the quark fields, on the other hand, cannot be directly put on a computer 
since quark fields are represented by anticommuting Grassmann numbers. 
Instead one uses the identity^] 


II d( ln dc ln 


exp 



l (U)q r) 



(27) 

(28) 


where D represents the lattice Dirac operator and <p n represents a bosonic 
field with 4 Dirac and 3 color indices to rewrite 


-^QCD — f dUnn f | d(f)^ l d(j) n exp ( *Sgluon ^nDnUU^m .(29) 

J n[i J n \ n,m / 

with which the fundamental variables U n ^ s and 0 n ’s are all bosonic. 

The inverse D^^iU) is a non-local quantity. A change of a (f) n spreads 
across the lattice through the inverse. Therefore preparing a trial configura¬ 
tion whose acceptance can be controlled is not straightforward. A number of 
methods were developed in the 80’s including the micro canonical plIHo] and 
Langevin j6F methods. The latter was also explored by the group at Cornell 
including Wilson [62]. The standard method has settled on the hybrid Monte 
Carlo (HMC) method proposed in 1987 [63]. which we now discuss. 

3 Strictly speaking, this identity requires positivity of the Hermitian part of ma¬ 

trix D. We shall not go into this technical detail and mention only that this can 

be guaranteed. 
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The first step of HMC is to introduce an SU(3)-algebra valued momentum 
P nil conjugate to J7 nAt , and rewrite the path integral (29) of full QCD as a 
partition function of a fictitious classical system of J7 n/ /s and P n ^s as 


^qcd = / n dPnn n dU nfJ , / , (30) 

n/i n/i n 

where H is a fictitious Hamiltonian defined by 

H =\Y j 1* - Sgluo„(f/) + £ ct>lD-} m (U)cf> m . (31) 

nfi n,m 


We now wish to generate a set of configurations distributed according to 
the weight exp(— H)/Zqcb by a Monte Carlo procedure. For this purpose, 
one introduces a fictitious time r conjugate to the Hamiltonian H. Starting 
with a given configuration of U n[x and <fi n at r = 0, we solve Hamilton’s 
equations for the canonical pair, 


— «7nA, 


d 


-p = 
dT ntl dU, 


Sgl uon (u) 


d 


+ E ^D-^—D^D-^U)^ 

n,k,l,m ° Un » 


(32) 


(33) 


over some interval of r. The configuration at a final time r is used as a new 
configuration in the Monte Carlo procedure. 

In numerical implementations, a continuous fictitious time evolution is 
discretized with a finite step size Sr. Since the Hamiltonian is no longer con¬ 
served, the configuration generated after a number of steps r n = nSr suffers 
from a bias. This is corrected by accepting or rejecting the generated con¬ 
figuration according to the Metropolis probability P acc = min(l, exp(— SH)) 
where 5H = H(P(T n ),U(r n )) — H(P(0),U(0)) is the difference of Hamilto¬ 
nian over the trajectory of length r n = nSr. 

Hybrid Monte Carlo is an elegant method. It is (i) exact, (ii) allows control 
of acceptance, i.e., the probability that a trial configuration is accepted, 
through the magnitude of Jr, and (iii) solving Hamltion’s equations can 
be executed in a parallel fashion. However, at every discretized step Sr in 
updating the fields according to the Hamilton’s equation, it requires the 
inverse of the lattice Dirac operator of form x n = fc> r given 

Unii and 0 n . 

The inverse can be obtained by solving the lattice Dirac equation, 


'Y^Dn,m{U)x m = 4>n- (34) 


This is a linear equation for a large but sparse matrix D n?m , which can 
be obtained by iterative algorithms such as the conjugate gradient method. 
The number of iterations N[ ter needed to reach an approximate solution is 
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controlled by the condition number k(D) of the matrix D. It is given by 
k(D) = A max /A m i n where A max and A m i n are the maximum and minimum 
value of the eigenvalues of D. Typically, the iteration number 7V iter is inversely 
proportional to the condition number. Since the minimum eigenvalue for the 
lattice Dirac operator is of the order of quark mass m q a in lattice units, and 
the maximum eigenvalue is 0(1), one finds N- lter oc l/m q a. 

Among the six types of quarks known experimentally, the lightest up 
and down quarks have masses of the order of a few MeV. This is three 
orders of magnitude smaller than the typical hadronic scale of 1 GeV. The 
condition number of the lattice Dirac operator for these quarks is large, and 
therefore, hybrid Monte Carlo simulation slows down considerably as one 
tries to approach the physical values of quark masses. 


3.5 Physical point calculation 


Lattice QCD simulations including quark effects through hybrid Monte Carlo 
algorithm developed rapidly toward the end of the 90’s. By the turn of the 
century, there was enough experience accumulated to empirically estimate 
how much computing is needed to calculate observables with a quotable 
error for a given quark mass. In Fig. [8] is shown a typical plot [64], taking 
the case of Nf = 2 flavor simulations for a lattice box of a size 3 fm, a 
minimum value to contain a hadron such as a nucleon within the box. The 
vertical axis shows the amount of computing in units of Tflopsxyear, i. e., 1 
year of running a computer executing 1 Tflop= 10 12 floating point operations 
per second. The horizontal axis is the ratio of i r meson to p meson masses, 
which varies with up and down quark masses and equals m^/rrip = 0.18 
in experiment. The three curves correspond to the inverse lattice spacing 
1/a = 1,2,3 GeV, with 1/a = 2 GeV or larger being needed for a reliable 
continuum extrapolation. We observe a sharp rise of the curves toward the 
physical value of m^/nrip = 0.18. This is a critical slowing down of the hybrid 
Monte Carlo algorithm, primarily arising from the slowdown of the Dirac 
solver toward m q a —>> 0. 

The rapid increase presented a major problem for lattice QCD simula¬ 
tions. Without overcoming this problem, one could only compute at quark 
masses much heavier than the physical values. The results had to be ex¬ 
trapolated to the physical point, but such an extrapolation involved large 
systematic errors because of the existence of logarithmic terms of the form 
^ m q log m q /A expected at vanishing quark mass. 

The difficulty was resolved in the middle of the 2000’s [65.. Since the force 
in the Hamilton’s equation (33) involves there are contributions 


coming from the short-distance neighborhood of the link np being updated 
and from those further away. It is possible to rewrite the quark determinant 
det D(U) such that these two types of contributions are separated. It was 
shown in [65] that, decomposing the lattice A into a set of sub-lattices A^i = 
1, 2, • • *, one can write 


det D{U) = ]3 det Di(U) ■ det R(U), 


(35) 
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Fig. 8 Computational cost for generating 1000 independent configurations for 
Nf = 2 flavor full QCD as a function of m n /m p whose physical value equals 
0.18. The three curves correspond to the lattice spacing a -1 = 1, 2, 3 GeV. Uni ts i s 
Tflopsx years, i.e., 1 year of time with the average speed of 1 Tflops. From Ref. [B3] . 


10 

1 


0.1 


1 2 3 4 d 5 

Fig. 9 Magnitude of the force term F^ on (k = 0), (k — 1), and F^(h = 2) 
as a function of distance from sub lattice boundary. From Ref. [65] . 



where Di(U ) is the Dirac operator restricted to a sub-lattice and R couples 
sites belonging to different sub-lattices. Rewriting each determinant factor 
on the right hand side as a bosonic integral, one can write the force term as 

Fn» = < u ° n + F™ + F®, (36) 

where F^ on is the term coming from the gluon action, those arising 

from Di(Uys, and hence represents short-distance contributions, and F^ is 
the term from R(U) giving long-distance contributions. 

In Fig. [9] we show the magnitude of the three terms observed in a hybrid 
Monte Carlo run [65]. There is a clear separation of magnitude for the gluon, 
and UV and IR quark contributions. Therefore, the step size Atir for the 
IR part F^ can be taken much larger than Arpv f° r the UV part, which in 
turn can be taken much larger than Sr g \ uon for the gluon part. An evolution 
with different step sizes for different parts of the force can be realized by a 
multi-time step integrator originally developed in [66 ] . 

Since inverting the long-distance part of the operator R is the compu¬ 
tationally most intensive, using different step sizes for the three parts of 
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the force lead to acceleration of the evolution by a large factor of order 
Srm/Sr g \ uon ^ 10 — 50. Simulations incorporating such an acceleration tech¬ 
nique were first carried out in the late 2000’s m, and reached the physical 
quark masses for the up and down quarks. 

The separation of UV and IR parts of the quark force can be realized in 
a different manner [68] by rewriting the Dirac operator for a quark mass mo 
as a product of ratios of successively heavier masses, 


N 

D(U, to 0 ) = ~[[(D(U,m i )D-\U,m i+1 )) ■ D(U,m N ), (37) 

z=0 


where rrq+i is taken larger than ?rq. Qualitatively speaking, one is separat¬ 
ing the contributions to the force according to the eigenvalues in the range 
[ra*, rai+i] with m^+i equal to the largest eigenvalue. 

A somewhat different method to accelerate the hybrid Monte Carlo al¬ 
gorithm is provided by the rational hybrid Monte Carlo algorithm [69]. One 
writes 


det D = 


det D 1/n 


/ n 

JJ d^dcj) 1 exp 
i= 1 


n 

i= 1 


(38) 


for some positive integer n, and applies a rational approximation to the 
fractional power x ~ x ! n . If k is the condition number of D, the n-th root 
D l / n will have a smaller condition number k}/ u . The magnitude of the force, 
compared with that of the hybrid Monte carlo method, will be reduced by 
a factor Hence the step size can be increased by the corresponding 

factor, leading to an acceleration of the molecular dynamics evolution. 

A final comment concerns the iterative inversion of the lattice Dirac oper¬ 
ator (34). Since the increase of the number of iterations toward small quark 
masses is the cause of the slowdown of the hybrid Monte Carlo algorithm, an 
ultimate optimization of the algorithm would be realized if one could remove 
the critical slowdown. This has recently been achieved by understanding the 
modes corresponding to the small eigenvalues of the lattice Dirac operator. 
One can either construct these modes explicitly and ’’deflate” (z.e., remove) 
them from the inversion CO], or employ multi-grid techniques to adaptively 
generate and incorporate these modes in the inversion E3MI. 


Over the years, the improvements as described here plus the development 
of an increasingly powerful computers have made it possible to carry out 
lattice QCD calculations at the physical quark masses. Such calculations 
are now routinely done. This is an impressive achievement. It also carries an 
aesthetic appeal; with the physical quark masses, one is no longer simulating, 
but rather calculating, the physical processes of the strong interactions as 
they are taking place in our universe. 
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4 Physics results 

4.1 Light Hadron spectrum 

4-1.1 Mass spectrum of light mesons and baryons 

Since the masses of hadrons are dynamical quantities, whether lattice QCD 
can quantitatively explain the experimentally known mass spectrum pro¬ 
vides a stringent test of the validity of QCD at low energies. Furthermore, 
the success of such a calculation forms the basis on which the reliability of 
lattice QCD predictions of other physical quantities are to be built. For these 
reasons, the calculation of the hadron mass spectrum, in particular those of 
the ground states of light mesons and baryons composed of up, down, and 
strange quarks, has been pursued since the beginning of lattice QCD and 
numerical simulations. 

A precise determination of the mass spectrum has to control a number of 
sources of errors. They are (i) statistical error due to the Monte Carlo nature 
of the calculation, (ii) systematic error due to extrapolation to the physical 
quark masses, (iii) systematic error due to finite lattice sizes, (iv) systematic 
error due to finite lattice spacings. In addition, until the late 90’s, most 
calculations were carried out with the so-called quenched approximation in 
which effects of quark-antiquark pair creation and annihilation are entirely 
ignored by dropping the quark determinant det D from the path integral. This 
was due to a high cost of including the quark effects into the simulations. 

The first systematic attempt at a precision calculation of the light hadron 
spectrum was carried out by Weingarten and his colleagues within the quenched 
approximation in 1993 m- The calculation took 1 year of dedicated com¬ 
puter time on GF11 computer developed by him at IBM. This was a landmark 
calculation in that a controlled extrapolation to the physical quark masses, 
to infinite volume and to zero lattice spacing, was systematically attempted. 
The masses of 7 r, p and K mesons were employed to fix degenerate up and 
down, and strange quark masses. The masses of K* and f> for the meson 
octet, nucleon N for the baryon octet, and A, A*, E* and Q for the baryon 
decuplet were obtained as predictions. They are plotted by filled circles in 
the left panel of Fig. [lO] The masses of A , £ and E of the baryon octet were 
not reported. The horizontal bars show the experimental values. We observe 
that the calculated values are consistent with experiment within one stan¬ 
dard deviation. However, for baryons, a sizable magnitude of the errors of 
0(10%) make it difficult to conclude if there is a precise agreement. 

A definitive calculation of the hadron spectrum within the quenched cal¬ 
culation followed in 2000 [73]. This work took half a year of CP-PACS com¬ 
puter which was 55 times faster than GF11. The results for the light meson 
and baryon ground states are shown in the right panel of Fig. [lOj As one 
clearly sees there, the quenched spectrum systematically deviates from the 
experimental spectrum. If one uses the K meson mass tuk as input to fix 
the strange quark mass (filled data points), the vector meson masses ijik* 
and m$ are smaller by 4% (4a) and 6% (5<r), the octet baryon masses are 
smaller by 6% - 9% (4 — 7a), and the decuplet mass splittings are smaller 
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Fig. 10 The left panel shows the res ult of the first systematic quenched light 
hadron mass spectrum calculation [73] published in 1993. The right panel shows 
the definitive quenched result m reported in 2000. 


by 30% on average. Alternatively, if the (j) meson mass is employed (open 
circles), uik* agrees with experiment within 0.8% (2a) and the discrepancies 
for baryon masses are much reduced. However, mx is larger by 11% (6cr). In 
other words, there is no way to match the entire spectrum beyond 5 to 10% 
precision in quenched QCD. 

The CP-PACS calculation heralded the end of the era of quenched cal¬ 
culations. Efforts toward full QCD simulations, which had already been tak¬ 
ing shape in the late 90’s, intensified. A first systematic calculation in full 
QCD was made by the CP-PACS Collaboration with dynamical up and down 
quarks (Nf = 2) [75] . With the algorithmic development in the middle of the 
2000’s, which we described in Sec. |3.5| calculations around the physical quark 
masses became possible. The PACS-CS Collaboration carried out Nf = 2 + 1 
flavor calculations in which the strange quark mass was taken close to the 
physical value and the pion mass was decreased down to m n = 155 MeV as 
compared to the physical value of 135 MeV m ■ 

Finally, the BMW Collaboration carried out an infinite volume and con¬ 
tinuum limit extrapolated calculation in 2008 m • We reproduce their results 
in Fig. 11 There is good agreement with the experimental spectrum within 
the errors of 2-5% except for A for which the error is 13%. 

Some of hadrons, including p and A, are resonances which decay for 
physical quark masses and infinite volume. A general relation between energy 
levels of two body states at finite volume and scattering phase shift, from 
which resonance parameters can be extracted, was established by Liischer 
some time ago m The BMW Collaboration assumed the Breit-Wigner form 
for the phase shift, and used this relation to correct the measured masses for 
the effects of decay coupling. 

A more rigorous approach in which the phase shift is first extracted from 
measured energy levels and Liischer’s formula, and resonance masses and 
widths are subsequently determined from the phase shift, was pioneered in 
[83] . Because of severe computational cost, it will take more time before 
physical point calculations can treat resonances numerically precisely in this 
way. 
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Fig. 11 Light hadron mass spectrum in Nf = 2 + 1 flavor QCD EH reported in 
2008. Masses of 7r, K and E are used to fix the degenerate up-down and strange 
quark masses as well as scale. Boxes of shaded grey show resonance masses and 
widths. Color online. 


4-1-2 Iso spin splittings of light hadron masses 

The isospin multiplets of hadrons exhibit small mass differences of a few 
MeV. These tiny effects are nonetheless very important to understand our 
universe. For example, neutron with a mass m n = 939.565379(21) MeV is 
heavier than proton of a mass m p = 938.272046(21) MeV by m n — m p = 
1.2933322(4) MeV. Because of this difference, neutron undergoes a /3 decay 
with a mean life of 880.3(1.1) sec, which has important consequences in nu¬ 
cleosynthesis after the Big Bang, and the composition of nuclei as we see 
them today. 

Since isospin mass splittings arise from the mass difference of up and down 
quarks, Nf = 1 + 1 + 1 simulations with different masses for up, down, and 
strange quarks become necessary. Once we take the quark mass difference 
into account, one also likes to consider the electromagnetic effects, since the 
magnitude of the effect, expected at the order of a few MeV, is similar to 
those arising from the up-down quark mass difference. 

A pioneering work on isospin breaking effects was carried out in the mid- 
90’s [78]. The RBC Collaboration rekindled interest by pointing out its im¬ 
portance m • Both quenched QED m and full QED [801181) have been 
attempted. The most recent calculation has been reported by the BMW 
Collaboration [82]. They carried out Nf = 1 + 1 + 1 + 1 QCD and QED 
simulations with independent up, down, strange, and charm quark masses 
and three values of QED coupling for a variety of lattice sizes and spacings. 
A careful analysis of finite size effects due to the infinite range of photon was 
made. The infinite volume and continuum limit extrapolation was carried 
out. The final result is 


m n -m p = +1.51(28) MeV (39) 

in good agreement with experiment. Treating the hadron mass differences to 
first order in m u — nid and QED coupling <a, they could separate QCD and 
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QED contributions with the result, 

m n — m v = +2.52(30)qcd — 1.00(16)qed MeV. (40) 

There is a delicate cancellation between the QCD and QED effects before 
the final number settles on the experimental value. 


4.2 Fundamental constants of the strong interaction 
4-2.1 Quark masses 

The determination of quark masses is a very important consequence of hadron 
mass spectrum calculations with lattice QCD. Since quarks are confined un¬ 
like leptons, lattice QCD provides the only reliable way for finding the values 
of this set of fundamental constants of our universe. 

In Table [4] we list representative lattice results as well as estimates by 
Particle Data Group over the years. Even as late as 1998, the Review of 
Particle Physics listed quite a wide band of values as shown in this table. 
This exemplifies how uncertain quark masses were only a decade and a half 
ago. Two years later, the CP-PACS quenched spectrum calculation narrowed 
the range considerably. However, the limitation of quenched QCD is manifest 
in a large discrepancy of the strange quark mass depending on the input. 

The Nf = 2 full QCD calculation [75] including dynamical up and down 
quarks but treating strange quark in the quenched approximation showed 
that (i) quark masses are significantly smaller than indicated by the quenched 


Table 4 Light quark masses in MeV units in the MS scheme at /a — 2 GeV. 


year 



action 

m u +m d 

2 

m u m d 

m s 

quenched QCD 




^ tyik input 

2 ) input 

2000 CP-PACS HU 

Wilson 

4.57(18) 

- 

115.6(2.3) 1} 

143.7(5.8) 2) 

Nf = 2 QCD 





^ ttik input 

2 ) m $ input 

2000 CP-PACS US] 

Wilson 

O ^r+0.14 

O.4tO_ 0- 20 

- 

89+g 

90t? o 2 > 

Nf - 2 + 1 QCD 

2010 MILC 
2010 BMW 
2012 RBC/ 
UKQC1 

84 

85 

D [ 

2Z] 

staggered 

Wilson 

domain 

wall 

3.19(18) 

3.469(67) 

3.37(12) 

1.96(14) 4.53(32) 

2.15(11) 4.79(14) 

89.0(4.8) 

95.5(1.9) 

92.3(2.3) 

Reviews of Particle Physics 

1998 PDG IBS] 
2012 PDG [HZ] 



2-6 
q k+0.9 
—0.2 

1.5-5 3-9 

2.3ig;£ 4.8 ±g;* 

60-170 

95 ±5 
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Table 5 Heavy quark masses in GeV units in the MS scheme at fi — mu- 


year 



action 

m c 

m h 

m t 

A/ = 2 + 1 QCD 

2010 HPQCD [89] 

staggered 

1.273(6) 

4.164(27) 

- 

1998 PDG 
2012 PDG 

|86J 

87 



1.1-1.4 
1.275(25) 

4.1-4.4 
4.18(3) 

173.8(5.2) 

173.07(89) 


results, and (ii) the discrepancy of strange quark mass depending on the input 
almost disappears. 

The recent results from Nf = 2 + 1 full QCD listed in Table [4] covers 3 
types of fermion actions. All three calculations carry out infinite volume and 
continuum extrapolations, albeit the degree of sophistication differs among 
the three. The separate values of up and down quark masses are estimated 
using additional input on isospin breaking such as the K° — K + mass dif¬ 
ference and estimations on electromagnetic effects. The three sets of results 
are reasonably consistent. The 2012 Review of Particle Physics values reflect 
this progress. 

In Table [5] we list the masses of charm and bottom quarks as determined 
by lattice QCD, together with those in Reviews of Particle Physics over the 
years. We also list the value for the top quark for completeness. 

Lattice QCD is rapidly moving into the era when all three light quarks are 
treated independently, and dynamical effects of charm quark is also included. 
It will be soon that such calculations yield a direct calculation of each of the 
quark masses with a few % error. 

4-2.2 Strong coupling constant 

The value of the strong coupling constant a s (/i) = g 2 (/a)/An defined in terms 
of the QCD coupling g{p) at some prescribed scale fi is another fundamental 
constant of our universe, on a par with the fine structure constant a = e 2 /47r 
for electromagnetism with the famous value a -1 = 137.035999074(44). 

In Fig [12] we plot the determinations obtained from scattering and decay 
data combined with perturbative QCD as of Fall 2013 [88] . and compare them 
with lattice QCD values from recent Nf = 2 + 1 calculations. By convention, 
the scale is taken to be the Z boson mass mz = 91.1876(21) GeV, and five 
flavors of quarks excluding top is incorporated in the running of the coupling. 
Lattice QCD determinations employ experimental quantities at low energies 
such as hadron masses and decay constants for fixing the scale. The methods 
employed range from step scaling function using the Schrodinger functional 
by PACS-CS [67] . current 2-point function and Wilson loops by HPQCD 189] . 
to static quark antiquark potential by Bazavov et al [HO]. The vertical band 
corresponds to the average over the four experimental determinations [88] . 

ai 5) (M z ) = 0.1183(12). 

The determinations with perturbative QCD still show some scatter and 
relatively large errors, and so do lattice QCD determinations depending on 
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Fig. 12 Strong coupling constant ai 5 \Mz) obtained from experiment and from 
lattice QCD. Here experiment means determinations from scattering and decay 
data combined with perturbative QCD. For further explanations, see text. 


the method, the ones from current 2-point function and Wilson loops having 
the smallest error. There is consistency with experiment at 1% level, and the 
precision of lattice determinations will steadily improve. 


4.3 CP violation in the Standard Model 


4.3.1 CKM matrix elements 


The connection between six types of quarks and CP violation is a salient fea¬ 
ture of the Standard Model. The complex phase characterizing CP violation 
is encoded in the Cabibbo-Kobayashi-Maskawa matrix V which appears in 
the weak interaction Lagrangian as 


Lint = ig (u,c,t) L 



wy + h.c. 

L 


In the Wolfenstein parametrization, the matrix takes the form, 


(41) 


f V ud ,v us ,v ub \ / l-X 2 /2,X,AX 3 (p-i V ) \ 

V= V cd ,V cs ,V cb = —A, 1 — A 2 /2, AX 2 +0( A 4 ), (42) 

\V td ,V ts ,V tb ) \- 1A 3 (1 p i:///).-,1A 2 . 1 ) 

with three real parameters A, A, p and an imaginary term irj characterizing 
CP violation. 

Constraining the CKM matrix requires matching experimental data on 
weak processes of hadrons to the theoretical expressions for the transition am¬ 
plitudes which follow from the interaction Lagrangian above. Since quarks in¬ 
teract strongly according to QCD, the weak transition amplitudes are dressed 
by QCD corrections. These corrections can be calculated by lattice QCD. 
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One of weak processes which plays an important role is the CP violating 
mixing of neutral K° and K mesons. The experimentally measured state 
mixing amplitude e can be written as 

e = CB k {r]iS(x c )Im(\ 2 c ) + r) 2 S(x t )Im( Aj) + 2r] 3 S(x c , x t )lm(X c X t )} .(43) 

Here C, 771,2,35 S(x,y), S(x) = S(x,x) are known constants and functions, 
x q = rriq/m^Y with m q and rriw the mass of quark q and W boson, and 
A q = Vq S V q d is the product of CKM matrix elements. The K meson bag 
parameter Bk is defined as 


B k 


(K°\(s^d) L (s^d) L \K°) 


8 

3 


Ik 171 


2 

K 


(44) 


where the expectation value is taken for QCD, and Bk is its renormalization 
group invariant value. We see then that the precision with which we can 
compute Bk directly reflects in the precision in the determination of the 
CKM matrix elements through the products X q = V* s V q d- 

Similarly, for the mixing of B q and B° q mesons with q = d or s, the 
oscillation frequency Am q is proportional to the B q meson decay constant 
fB q and the bag parameter Bs q defined by 




Vk m k 


(45) 


In addition, the decays of B mesons in the leptonic (B tv) and semi- 

leptonic (such as B niv and B —>> D*£v) channels are used to constrain 
the matrix elements | | and | |. 

The constraint on the CKM matrix is usually written on the complex 
(p, rj) plane defined by 


P + iV = ( 1 - y ) (P + irj) ■ 


(46) 


In Fig. 


13 we show the latest result compiled by the UTfit Group [91 . The 

is similar. The inside of various regions 


result by the CKMfitter group [92] 
are allowed from each of the constraints. There has to be a common overlap 
region if the Standard Model is to be consistent with experiment. Fitting 
data to the Standard Model, UTfit finds 


p = 0.132(23), 77 = 0.351(13). (47) 

For each region in Fig. [13] the width of the allowed band represents ex¬ 
perimental uncertainties as well as those of lattice QCD determinations. In 
the two left columns of Table [6] we list the experimental inputs, and their 
percentage errors, used for constraining the allowed regions. The quantities 
and percentage errors listed in the two right columns are the QCD matrix 
elements which are needed to convert the experimental inputs to constraints 
in the (p, rj) plane. 






33 



Fig. 13 Latest c ons traints on the CKM matrix element expressed in the (p, 77 ) 
plane. From Ref. [92] . Color online. 


The experimental values are taken from Reviews of Particle Physics (9] 
and a compilation of heavy flavor data by the Heavy Flavor Averaging 
Group [93] , The lattice QCD values for hadronic matrix elements are taken 
from a recent compilation by the Flavor Lattice Averaging Group (FLAG) [94] . 
We select the values quoted for Nf = 2 + 1 calculations that passes FLAG 
criteria for the control of systematics errors including the continuum extrap¬ 
olation. In some cases, only a few calculations are available, in which cases 
we quote the original references. 

We observe that B meson related quantities still has a significant mar¬ 
gin for improvement, both on the experimental and lattice QCD sides. The 
SuperKEKB experiment, which will start in a few years, is expected to re¬ 
duce the experimental error by another order of magnitude. Improvements in 
lattice QCD values should keep pace. We should also remark that the large 


Table 6 CP violating observables and relevant matrix elements. Percentage values 
of errors are also listed. See text for explanations. 


observable 

PDG ED/ 

HFAG [93] 


matrix 

element 

FLAG M 


e 

2.228(11) x 10 -3 

0.49% 

B k 

0.7661(99) 

1.3% 

Am s /Arrid 

34.8(2) 

0.60% 

f B s \/ b B s 

f B d s/BT d 

1.268(63) [SU 

5.0% 

B[B tv) 

1.14(27) x 10 “ 4 

24% 

f B ( MeV) 

190.5(4.2) 

2 . 2 % 

B[B 

0.38(2) x 10 “ 4 

5% 

ac b ^ 

2.16(50) |SB||57) 

23% 

B -> D*£v 

35.85(45) x 10 “ 3 

1.3% 

F b ^ d ’{ 1 ) 

0.906(13) [98| 

1.4% 


integral over 16GeV 2 < q 2 < g^ ax 
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error in the B decay measurement affects the determination of V^>, which in 
turn broadens the band for e. The determination of Bk itself is already quite 
precise. 


4-3.2 CP violation in the two pion decays of K meson 


Historically, CP violation was discovered through an observation of K mesons 
decaying into two pions in 1964. The strength of direct CP violation relative 
to that in the state mixing is measured by the ratio 


e f uj [ Irm4.2 Irm4o 
e V%\e\ ReAo 


(48) 


where Ai denotes the decay amplitude with isospin I in the final state, and 


i ReHo 
R eA 2 


(49) 


Two heroic experiments, NA48 at CERN 99 and KTEV at FNAL |100| . 
spanning two decades from the 80’s to early 2000’s, measured e'/e with the 
result, 


e 


/ 


e 


18.5(7.3) x 10- 4 , NA48 
20.7(2.8) x 10- 4 , KTeV 


(50) 


It has also been a long standing puzzle that the amplitude for the 1 = 0 final 
state is sizably larger than that for I = 2, namely cj - 1 ~ 22. This is called the 
AI = 1/2 rule. Whether the Standard Model can successfully explain these 
features of K —>> 2tt decay has been a major problem in particle physics. Since 
the issue boils down to calculating the strong interaction corrections to the 
effective weak interaction Hamiltonian, this has been an important challenge 
to lattice QCD since the 80’s. 

There are three obstacles to a successful calculation of the K —)> tttt am¬ 
plitudes in lattice QCD. The first obstacle is chiral symmetry. Chirality plays 
an essential role in weak interactions. The effective 4-quark weak interaction 
Hamiltonian obtained at a lattice cutoff scale of a few GeV starting from a 
much higher weak interaction scale has a definite chiral structure [101] . Thus 
one has to employ a lattice fermion formulation which has chiral symmetry. If, 
on the other hand, one uses non-chiral formulations such as Wilson’s, one has 
to carefully control chiral symmetry violation effects under renormalization. 
The former option was not available until the late 90’s when domain-wall 
and overlap formulations were proposed. The latter problem was successfully 
resolved for the Wilson fermion action only recently, with an enticing con¬ 
clusion that the renormalization structure is the same as in the continuum 
for K 7T7T decay [10211103| . 

The second obstacle stems from the fact that, in the Euclidean Green’s 
function for the K tttt transition by the weak Hamiltonian, the two-pion 
state with zero relative momentum, being the state with lowest energy in this 
channel, dominates for large times. This contradicts the physical kinematics 
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of the decay; the two pions decaying f rom a K mes on at rest should have an 
equal and opposite momentum p = \/m 2 K /4 — m\. Thus a naive calculation 
does not yield physical results. A resolution was found in 2001 |104j . One 
makes a calculation for a finite lattice volume chosen such that the energy 
of the two pion state matches the energy of the K meson. The physical 
amplitude for infinite volume can then be obtained by the following formula, 


I \hys(K 7T7r)| 2 = 87 [ P ^~ + 4^^} I I- H W M, ot | 2 .(51) 


Here p is the pion momentum, S n7r (p) the elastic tttt phase shift at momentum 
p, q = pL/(2ir), and <p(q ) is defined by 



(52) 


The third issue is specific to the I = 0 channel and computational in char¬ 
acter. In this channel, there are diagrams with disconnected quark loops, e.g., 
quarks from the pions annihilate themselves. In addition, the so-called Pen¬ 
guin diagrams in which a pair of quarks from the weak Hamiltonian forms 
a loop are also present. These diagrams suffer from large statistical fluctu¬ 
ations, rendering the statistical average difficult. This problem is gradually 
being overcome with efficient algorithms for computing disconnected and 
Penguin contributions, and with increase of computing power which allows 
a large number of Monte Carlo ensembles. 

Recently there has been significant progress assembling these develop¬ 
ments together. The RBC Collaboration has developed applications of the 
domain-wall formulation of QCD having chiral symmetry. It has succeeded 
in calculating the 1 = 2 amplitude for the physical pion mass [ 1051 . Their 
result, obtained at a lattice spacing of a ~ 0.14 fm, is given by 


Re^4 2 = +1.381(46) st at(258) syst x l(T 8 GeV, 
Im A 2 = -6.54(46) stat (120) syst x l(T 13 GeV. 


(53) 

(54) 


The real part is in good agreement with experiment: ReA. 2 Xp = 1.479(4) x 
10“ 8 GeV. 


The RBC Collaboration has been attacking the much more difficult 1 = 0 
channel, using a special G-periodic boundary condition to force the pions 
to carry momentum so that an energy matching is realized between the K 
meson and 2 pions. Preliminary results have been presented at Lattice 2014 
Conference this year. 

Another group at Tsukuba has been pursuing the problem using the Wil¬ 
son fermion formulation m- The renormalization structure for the operator 
relevant for K tttt decay turned out to be the same as in the continuum 
except for the mixing with dimension 3 operator 575 d, which is subtracted 
non-perturbatively. Their calculation so far does not achieve physical kine¬ 
matics. The 7 r meson mass is artificially taken large so that 2m^ = mx is 
satisfied. Nonetheless they reported an encouraging first result for I = 0 
channel as well as for / = 2 at Lattice 2014 Conference. 
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It is hoped that progress in the near future brings definitive results on 
the 1 = 0 amplitude Aq. This will put an independent constraint on the CP 
violating phase fj in Fig. 13 as the ratio e'/e is proportional to fj. 


4.4 Quark-gluon matter at high temperature and density 

Confinement of quarks and spontaneous breakdown of chiral symmetry are 
both dynamical consequences of QCD. A very interesting question then is 
how these properties may or may not change if parameters (or dials) external 
to QCD are varied. One of important dials is temperature, which increases 
toward the Big Bang in the Early Universe. Another dial is baryon number 
density, which has a large value in extreme conditions such as in the core of 
neutron stars. 

In both cases one is interested in an aggregate of hadrons, rather than 
individual hadrons, and wants to understand how its properties change at 
extremely high temperature or density. As we discuss below, a general con¬ 
clusion from lattice studies is that a gas of hadrons turns into a different 
state in which quark and gluon degrees of freedom becomes manifest. The 
novel state of matter is often called quark-gluon plasma (QGP). 

4-4-1 Phase diagram of QCD at finite temperature:analytical considerations 

The Euclidean formulation adopted for lattice QCD is well suited for studies 
of its properties at finite temperatures. If one considers a lattice with N t sites 
in the temporal direction, and imposes periodic and antiperiodic boundary 
condition for gluon and quark fields, respectively, the lattice path integral 
([5]) is equal to the canonical partition function 

Zqcd = Tr [exp (—Hqcd/T)] , (55) 

at a physical temperature 


T = 


N t a 


(56) 


This connection shows that methods developed for zero temperature stud¬ 
ies, including Monte Carlo calculations, are readily applicable to the finite 
temperature case. 

In order to discuss what happens as the temperature is raised from zero, it 
is helpful to consider a two-dimensional plane of the average up-down quark 
mass m u d and the strange quark mass m s , often dubbed Columbia Plot, as 
depicted in Fig. 14 |l06j . One can examine various limiting cases as follows. 


(i) Pure gluon theory: The top right corner in Fig. [l4] for m u d = m s = 00 
corresponds to the pure gluon theory without quark degrees of freedom. This 
case is important, nonetheless, since confinement is a dynamical consequence 
of the gluon fields. 
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N f = 2QCD SU(3) YM 



Np= 2 QCD SU(3) YM 



Fig. 14 Phase digram in 2+1 flavor QCD as a function of the degenerate u and d 
quark mass m u d and the s quark mass mn s . Left panel shows the case for a second 
order chiral transition for two-flavor QCD. Right panel shows the case for a first 
order two-flavor chiral transition. From Ref. [107| . Color online. 


The pure gluon theory possesses a center Z( 3) symmetry defined by the 
transformation 


UnA C^n4> C C ^(3), (57) 

for the sites n on some fixed time slice t. The corresponding order parameter 
is the Wilson loop winding around the space-time in the time direction at a 
fixed spatial site n defined by 


Q{ n) = Tr 


N t 


^(n,n t )4 


\n t = l 


(58) 


This operator, often called Polyakov loop, transforms under the center Z(3) 
symmetry as 


i?(n) -+ £J?(n). (59) 

Polyakov [108] and Susskind argued that the center Z(3) symmetry is 
intact at low temperatures with (i?(n)) = 0, but becomes spontaneously 
broken beyond a certain temperature T — T c with a non-zero expectation 
value (i?(n)) ^ 0. 

The connection with deconfinement is most easily understood in the fol¬ 
lowing way. The 2-point correlation function of Polyakov loops defined by 
+( n )+ m)) is connected with the free energy F(|n — m|) of a static quark 
at site n and an antiquark at site m by 

exp (— F(r)/T) = (i?(n)i?^(m)) , r = |n — m|. (60) 

Assuming the presence of a mass gap /+ the 2-point function on the right 
hand side is expected to behave for large r as 

(^(n^^m)) -+ (i?(n)) • (i?^(m)) + O (exp (— fir)) , r -+ oo. (61) 
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Depending on the expectation value (i?), one finds for r -> oc that 

( err , (Q) =0, T <T C 

F(r ) -» ^ (62) 

[ c +0(exp(-/x D r)), (j?) ^ 0, T > T c 

where a = T /a for T < T c , and c is a constant and = /i is a color electric 
Debye screening mass for T > T c . 

Analytical considerations indicate that this deconfinement phase transi¬ 
tion should be of first order m- If this is the case, the first-order transition 
will extend into the region of large but finite quark masses Em], terminating 
at a line of second-order transition as depicted in Fig. [l4j 


14 


(ii) Nf = 2 and Nf = 3 massless QCD: The top left corner of Fig. 
corresponds to Nf = 2 massless QCD with m U( i = 0 but m s = oo, and the 
bottom left corner to Nf = 3 massless QCD with m U( i — 0 and m s - 0. 

For Nf flavors of massless quarks, QCD is invariant under a global SU (Nf)® 
SU(Nf) chiral symmetry, which is spontaneously broken down to vector 
SU(Nf) symmetry with a non-zero vacuum expectation value of the order 
parameter N = (q n q n ) 7 ^ 0- The Nf — 1-plet of pseudo scalar mesons are the 
corresponding Nambu-Goldstone bosons. 

When one raises the temperature, thermal fluctuations tend to destabilize 
the chiral order parameter N. Thus one expects a phase transition restoring 
chiral symmetry at some temperature. A more detailed examination is pos¬ 
sible using an effective non-linear sigma model for the order parameter field 
= QnQn with h 3 — 1, • • • ,7V/ and renormalization group methods. The 
result Ena indicates that (i) for Nf = 3 or larger, chiral symmetry is re¬ 
stored through a first order phase transition, while (ii) for Nf = 2 the phase 
transition is of first or of second order depending on whether flavor singlet 
U( 1) axial symmetry is effectively restored or not at T > T c . 


(Hi) Connecting Nf = 2 and Nf = 3 massless QCD: One can interpolate 
between the Nf = 2 and 3 cases by changing the strange quark mass m s 
from oc to 0. If the Nf = 2 transition is of second order, the first order 
transition for Nf = 3 at m s = 0 has to change to a second order transition 
at a tricritical point at some m s = m c s . 


(iv) Nf = 3 symmetric QCD: In general, a first order phase transition is 
stable under symmetry breaking perturbations up to a critical value where 
it terminates with a second order phase transition. For the second order 
case, the phase transition immediately disappears if symmetry breaking per¬ 
turbations are turned on. Thus one expects a sheet of first order transition 
extending from the bottom left corner corresponding to Nf = 3 massless 
QCD to the interior of the phase diagram. The sheet should terminate at 
a line of second order phase transitions, which is expected to belong to the 
Ising universality class m- 

If the Nf = 2 transition at m u d = 0 is of second order, this line of second 
order transitions will hit the vertical axis at m s = m c a with a p ower law 


m ud oc (m c s - m s ) 5/2 mu, as depicted in the left panel of Fig. Il4l If, on the 
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other hand, the Nf = 2 transition is of first order, the line will go up to the 
top horizontal line (the right panel of Fig. 14). 


4-4-% Monte Carlo study of the finite-temperature phase diagram 

Lattice QCD simulations are carried out for a fixed temporal lattice size N t . 
One then regulates the temperature indirectly by varying the bare gauge 
coupling constant g since the continuum limit a = 0 is located at g g = 0, 
weaker couplings correspond to smaller lattice spacing, and hence to higher 
temperatures, and vice versa for stronger couplings and to lower tempera¬ 
tures. 

The basic tool for studying the phase diagram is finite size scaling theory 
developed in the late 60’s and 70’s ESI- This theory helps to analyze how 
singularities marking phase transitions develop from numerical data obtained 
at finite volumes. 


(i) Pure gluon theory: The first instance where finite size scaling method 
was crucially effective was the deconfinment transition of the pure gluon 
theory. Early simulations quickly found that the Polyakov loop expectation 
value exhibited a rapid increase from (i?) ~ 0 to non-zero values over a 
narrow range of temperature as expected. A subtle issue is if a rapid increase 
actually signifies a phase transition, and if so, whether it is a first order phase 
transition with a discontinuous jump of (f2) at T c or a second order phase 
transition with a continuous (i?) but having a singular derivative. 

The susceptibility of Polyakov loop is defined by 

Xo = JdJ2 (f?(n)f? f (m)) (63) 

n,m 

with L the linear extent of the system, and d the space dimension. For a 
finite volume, the susceptibility exhibits a peak. The infinite behavior of the 
peak height X^ ax distinguishes the order of the phase transition. According 
to finite size scaling theory, if one parametrizes the volume dependence by a 
power law, 


X £ ax (L) oc L a , oo, (64) 

the value a = d, i.e., the dimensionality of space, signals a first order transi¬ 
tion, while values a < d characterizes a second order transition in which case 
a is related to the critical exponents. Similar criterion holds for the volume 
dependence of the width of the peak. 

These criteria were utilized to establish the first order nature of the de¬ 
confinement transition in the pure gluon theory [ f 16] . in agreement with the 
analytical considerations. 
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(ii) Chiral transition with quarks: The chiral phase transition in the pres¬ 
ence of quarks has also be studied extensively. Ideally one would like to use 
a fermion formulation preserving chiral symmetry such as the domain-wall 
or overlap. They became available only in the early 2000’s, and are compu¬ 
tationally very costly, however. For these reasons, most of the calculations 
have utilized, and still use, the staggered fermion formulation and, to a lesser 
extent, the Wilson formulation. 

Broadly speaking, results accumulated to date are as follows. For Nf = 3 
degenerate flavors, the phase transition is consistent with being of first order 
for small quark masses. Increasing the quark mass, the transition weakens 
and terminates at a second order transition whose exponents are consistent 
with the Ising universality class (see Refs. (11711118] for early representative 
work and references). These calculations are made at most for N t = 6 lattices, 
and the continuum limit is yet to be taken 

For Nf = 2, the situation is more complicated. Early simulations were 
consistent with a second order phase transition both with the staggered ( 1201 
and with the Wilson P2B fermion formulations. For the staggered case, how¬ 
ever, the critical exponents did not come out consistent with either 0(4) nor 
0(2) values. There have also been simulations suggesting consistency with a 
first order transition |122| . 

As already pointed out in Ref. na, the order of the Nf = 2 chiral tran¬ 
sition is connected with the Ua{ 1) anomaly. Recently, a theoretical argument 
has been put forward that the anomaly effects disappear for certain sets 
of correlation functions if chiral symmetry is restored. At the moment, the 
order of the Nf = 2 transition is not settled completely. 


4-4-3 Thermodynamics with physical quark masses 

Physically, the crucial issue is where the point with physical quark masses lies 
on the phase diagram in Fig. [M] The staggered results are unanimous that 
it lies beyond the line of critical end points. Hence there is only a continuous 
crossover and no phase transition with a singular behavior. The basis for 
this conclusion includes an extensive study at physical quark masses with 
infinite volume and continuum extrapolations using susceptibilities of various 
physical observable |126| . 

Since the transition is a continuous crossover, the transition tempera¬ 
ture is not uniquely determined but depends on the quantity used. For ex¬ 
ample [127] . one finds T c = 147(2)(3) MeV if one uses the susceptibility 
peak of chiral order parameter, and T c = 157(4) (5) MeV from the inflec¬ 
tion point of the energy density. An independent calculation [129] reported 
T c = 154(9) MeV from 0(4) scaling analysis of chiral susceptibility. The 
results are consistent, and altogether indicate T « 150 — 160 MeV as the 
temperature range across which the physical chiral transition takes place. 

In the left panel of Fig. [l5] we show the energy density e/T 4 and in¬ 
teraction measure J/T 4 = (e — 3p)/T 4 , with p the pressure, in units of T 4 
calculated in the staggered quark formalism by two groups [T281H30] . They 
are obtained at the physical quark masses, and infinite volume and con¬ 
tinuum extrapolations are made. We observe very good agreement up to 
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Fig. 15 Thermodynamic quantities in the continuum limit in 2+1 flavor QCD as a 
function of temperature. Left panel shows the energy density e/T 4 and interaction 
measure I /T 4 = (e — 3p)/T 4 in units of T 4 , and right panel shows energy density 
e in units of GeV/fm 3 . From Refs. [128U130| . Color online. 


about 200 MeV beyond which a one standard deviation difference appears. 
Physically important is the feature that the Stefan-Boltzmann value for free 
gluons and quarks, esB /T 4 = 7t 2 (8 + 217V//4)/15 = 15.62 • • • for Nf = 3, is 
reached only slowly, with significant deviations remaining at T/T c ~ 2 — 3. 
This indicates that the quark-gluon matter is strongly interacting at these 
temperature ranges. The right panel shows the energy density in units of 
GeV/fm 3 . 

Experimental effort toward detection of quark-gluon plasma through heavy 
ion collisions in accelerators has been going on for a long time. It started with 
the Bevalac at Berkeley in the 70’s, the AGS at BNL and the SPS at CERN 
in the 80’s, with RHIC at BNL since 2000 and with the LHC at CERN most 
recently. 

Let us see how lattice QCD predictions compare with heavy ion ex¬ 
periments at RHIC and the LHC. Phenomenological estimates on the en¬ 
ergy density reached at the initial stage of the collision range from e ~ 
5.6 GeV/fm 3 for Au-Au collisions at RHIC with s/snn = 200 GeV [131 j, to 
e « 15 GeV/fm 3 for Pb-Pb collisions at the LHC with y / +Viv = 2.76 TeV [132] . 
Looking at the right panel of Fig. [15] we read out a temperature of order 
T ~ 200 MeV and 300 MeV, respectively, for these energy densities which 
are high enough for the collision product to be in the high temperature phase. 
For comparison, an experimental estimate of the temperature using the low 
Pt excess of direct photons that are supposed to come from the initial ther- 
malized state gave T = 221(27) MeV at RHIC [133] and T = 304(51) MeV 
at the LHC [134] . 

The initial fireball rapidly cools as it expands, and hadrons are formed 
once the temperature falls below the transition temperature T c . It has been 
known that the yield of various hadrons from the collision could be well 
fitted with statistical thermal distribution parametrized by a chemical freeze 
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Fig. 16 Schematic phase diagram of QCD on the temperature baryon density 
plane.. 

out temperature Tf. This temperature increases with the collision energy 
and saturates at Tj ^ 160 MeV at RHIC energies [135] , though slightly 
decreasing at LHC energies. 

Furthermore, the azimuthal anisotropy in the transverse hadron yields, 
quantified in terms of the elliptic flow v 2 and higher moments where v n = 
(cos ri0), is well described by hydro dynamical models with no or small vis¬ 
cosity and using the equation of state from lattice QCD. 

Thus, while small viscosity was a surprise, the experimental data are 
broadly consistent with the production of strongly interacting high temper¬ 
ature quark-gluon matter with the properties as predicted by lattice QCD. 
It is interesting to note that small shear viscosity close to the quantum limit 
p/s = 1/47r |136] is observed by lattice QCD calculations of transport coef¬ 
ficients in the pure gluon case [13711138] . 

From a heavy-ion collision point of view, the transition temperature and 
equation of state are only indirectly reflected in the characteristics of the final 
hadronic state. The moments of conserved charges such as electric charge Q, 
strangeness S', and baryon number B provide interesting quantities which are 
calculable in lattice QCD and are directly observable in experiments. Thus 
they have been attracting much interest recently [1391I140U141U142U143] . 

4-4-4 Dynamics at finite baryon number density 

Theoretical expectations on the phases of QCD at finite temperature T and 
finite baryon number density pB is depicted in Fig. [16] The chiral crossover 
transition at T c « 150 — 160 MeV continues into the finite pB region, and hits 
a second order critical point at some (Te,Pe) beyond which the transition 
turns into a first order phase transition. For sufficiently large baryon number 
density, one expects novel phases such as a color superconductor. 

Finite baryon number density can be introduced by adding a quark chem¬ 
ical potential p to the Hamiltonian. On a lattice, this is achieved by multi¬ 
plying the positive temporal hopping term of the quark action by exp (pa), 
and the negative hopping term by exp (—pa). 
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Solid quantitative results are still not available on the phase diagram. The 
main reason is that for non-zero quark chemical potential, the quark Dirac 
determinant det D(U,/jl) is complex, so that the Monte Carlo methods based 
on a probability interpretation of the weight det D • exp(— Sg(U)) no longer 
work in general. 

One can see the difficulty by defining the phase of the determinant by 


0(u, n) = -i log 


det D(U, n) 


| det D(U,./j)\’ 
and rewriting expectation values in the following way: 

(O exp (i6(U, fi)))n 


( 0 ) = 


(exp (i0(U, n))}. 


(65) 


( 66 ) 


where (-)|| means the average with respect to the phase quenched determinant 
| det D(U, fi)\. The quenched average of the phase factor in the denominator 
is a ratio of the two partition functions: 


(exp (id(U, n))) 


Z \\M 


W 

/ n dU niJ ,| det D(U, fi)\ exp(S g i uon {U)). 


(67) 

( 68 ) 


If the quenched average defines a statistical system with a free energy density 
/||(/i), one can write 


/ / -am ^ ( n/(M)-/||(M)\ 

(exp (i9(U, m)))|| = exp I-— — u - J , (69) 

with f(fi) the free energy density of the original system. The right hand side 
vanishes exponentially fast for large spatial volumes V —>• oo. This is one way 
of explaining the sign problem. 

A variety of methods have been devised and explored to overcome this 
problem. Besides the phase quenched calculation described above, they in¬ 
clude reweighting of the determinant |144j , analytic continuation from imag¬ 
inary chemical potential |145ill46j . Taylor expansion in powers of fi/T |147l 
1148] , canonical ensemble simulation [ 149 , complex Langevin simulation (150 , 
H5B, and others. These methods have yielded some success for not too large 
values of ji/T. The main problem, however, has been the difficulty of con¬ 
trolling errors in the results. The precise location of the critical point E and 
physics of finite density QCD for larger baryon number density is still largely 
open at present. 
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5 Conclusions 

Lattice QCD is a major contribution of Kenneth Wilson to physics. In our 
view, it is on a par in significance with his renormalization group theory 
of critical phenomena which won him Nobel Prize in Physics in 1982. Con¬ 
ceptually, it clarified how quantum fluctuations of gauge fields give rise to a 
confining force which is essentially distinct in its origin from the forces arising 
from exchange of a particle such as the electromagnetic force. At the same 
time, coupled with supercomputers, it opened a way to calculate the physi¬ 
cal predictions from its first principles, making possible detailed comparisons 
with experiment. 

In 2004, Kenneth Wilson delivered a talk entitled “The Origins of Lattice 
Gauge Theory” at International Symposium on Lattice Field Theory held at 
Fermi National Accelerator Laboratory |[T. Looking back on the development 
of lattice QCD, he said that “The lattice gauge theory was a discovery waiting 
to happen, once asymptotic freedom was established. ”, and went on to describe 
works of his contemporaries who could have preceded or shared the discovery. 
Nonetheless it is clear that, because of his previous studies leading to his 
encounter with statistical mechanics and phase transitions, he was in a unique 
position to be the first to grasp the deep significance of strongly coupled gauge 
dynamics in relation to confinement. 

In the same article, he commended on the vast progress in lattice QCD in 
the thirty years since 1974, the year of his seminal paper pQ, due “in part to 
improved algorithms, in part to increased computing power, and in part to the 
increased scale of the research effort underway today ”, but urged on that “this 
does not mean that the present state of lattice gauge computations is fully 
satisfactory. The knowledge and further advances that will likely accumulate 
over the next thirty years should be just as profound”. 

As if corroborating his view, in just a decade since then, the physical 
point computation was realized, and many beautiful physics results ensued, 
as summarized in Sec. [U 

The past decade, then, was a turning point of lattice QCD in our view; 
prior to this event, despite the premise that it provided a first-principle cal¬ 
culation, it remained a method of uncertain reliability, requiring chiral ex¬ 
trapolations which were difficult to control. We are now able to calculate and 
understand many properties of hadron states including masses and matrix 
elements, at the physical quark masses and controlling errors of calculations. 
The precision of the calculation has now reached the level of a few % or better 
in many quantities. Consequently, important constraints are now available on 
the CKM matrix elements and CP violation in the Standard Model. 

This is not to say that progress has been uniform in all fronts. In thermo¬ 
dynamics of QCD, while some basic quantities such as the transition tem¬ 
perature and equation of state have been calculated, the region of finite 
baryon number density is still largely unexplored. A major methodological 
breakthrough will be needed to extend our understanding to the entire phase 
diagram of QCD. 

One area which did not even exist at the time of Wilson’s talk in 2004 
is nuclear physics based on lattice QCD. In fact, while there were pioneer- 
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ing studies on H di-baryon in the late 80’s fl52H153| and nucleon scattering 
lengths in 1994 |154j . it is in 2007 that the two nucleon potential was first 
extracted from lattice QCD |155j . and in 2009 that the binding energy of 
Helium was calculated directly from Helium correlation function [156] . This 
is a challenging area in terms of physics as well as in calculational tech¬ 
niques; one has to deal with a small energy scale of 0(10) MeV, the number 
of Wick contractions for quark fields increases factorially fast with the mass 
number, and so does the statistical error of nuclei propagators for large time 
separations. 

In theoretical physics, theory and computing go hand in hand. Calcula¬ 
tions are indispensable in order to confirm the validity of theory, and even 
more so to explore the consequences of theory which help us better under¬ 
stand why our world works the way it does. Lattice QCD is a prime exam¬ 
ple of such a relationship between theory and computing. Kenneth Wilson 
clearly foresaw the importance and future potential of supercomputing in 
this connection, and many of his thoughts and vision in the early 80’s 
came to be realized since then. Looking toward future, as he prophesied [2], 
our understanding of the strong interactions will become more profound as 
the computing power increases toward exascale and possibly beyond in the 
decades to come. 
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